# Formalism to deal with discontinuous potentials in classical mechanics (hard wall, hard spheres)

It seems to me that Hamiltonian formalism does not suit well for problems involving instantaneous change of momentum, like particle collisions with hard wall or hard sphere gas model. At least I could not apply it straightforwardly to the simplest possible problem of 1D particle hitting a wall:

            │/ wall
│/
particle  │/
───o────────┼────────────> x
│/
│/
│/


My attempt was quite direct. I took the Hamiltonian to be

$$H = \frac{p^2}{2m} + U(x)$$

with potential $U$ defined as

$$U(x) = \cases{ 0, \; x < 0, \\[.5em] K, \; x > 0, \\[.5em] E, \; x= 0.}$$

where $E$ is the particle energy and $K > E$. Hamiltonian equations should read as

$$\cases{\dot x = \frac{p}{m}, \\[.6em] \dot p = - \frac{d U}{d x}.}$$

It is not hard to integrate the first equation, but my attempts to integrate the second one did not lead to any meaningful result (that's why I do not share them here, it was a complete failure).

So I ask whether it is possible to obtain the solution to the problem by directly integrating Hamiltonian equations in the form above, without relying on general mechanical theorems/principles like energy conservation? Or is such an approach completely unsuitable for the task?

If so, what is the general (and elegant) approach to such systems?

There exist a related question on PSE "Hamiltonian function for classical hard-sphere elastic collision", but the setting is more cumbersome.

• One way to go about it is to switch to using $x$ as the independent variable and $t$, i.e. the time of flight to get to a specific $x$, as a dependent variable. Then your discontinuous-in-x Hamiltonian becomes a step function kick.
– webb
Commented Mar 27, 2014 at 16:24
• If noone comes with an elegant framework which would avoid potential regularisation in the remaining bounty time, I will accept Qmechanic's answer. Commented Mar 31, 2014 at 16:19

## 5 Answers

Perhaps the simplest and most intuitive approach is to regularize the hard wall potential

$$V_0(x)~=~\left\{ \begin{array}{rcl} 0 &\text{for}& x<0 \cr\cr \infty &\text{for}& x>0\end{array}\right.$$

as

$$\lim_{\varepsilon \to 0^+} V_{\varepsilon}(x) ~=~V_0(x).$$

For instance, one could choose the regularized potential as

$$V_{\varepsilon}(x)~=~\frac{x}{\varepsilon}\theta(x).$$

This corresponds to constant velocity for $x<0$ and constant acceleration for $x>0$. Next write down a continuous solution for the position $x_{\varepsilon}(t)$ as a function of time $t$, say, for given pertinent initial conditions.

Finally, at the end of the calculation, one should remove the regularization $\varepsilon \to 0^+$ again, and check if the limit $$\lim_{\varepsilon \to 0^+} x_{\varepsilon}(t)$$ makes physical sense.

• Yes, this regularized potential works as a charm. It reverses the initial velocity and the penetration depth in the realm $x > 0$ is $\varepsilon E$, while the time spent here is $1/2 \, \pi \sqrt{\varepsilon m}$. Both these expressions go to zero with $\varepsilon \to 0^+$. I have one question though, how did you come up with this potential? Of course it is possible and probable to deduce it, but I bet you just had the right analogy from you experience. Is it so? Commented Mar 31, 2014 at 16:15
• @Yrogirg: Right, when dealing with idealized collisions, one often encounters infinite potentials and infinite forces. Of course, in reality, the colliding objects deform, and the forces involved are only finite in value. The regularization $\varepsilon$ can be viewed as imitating such (elastic) deformation. Commented Mar 31, 2014 at 16:30
• oops, I have awfully miscalculated the time. If somebody cares, the correct reasoning should be following: The distance traveled is $\varepsilon E_0 = v_0 (\tau / 2) - 1/2 \, \varepsilon (\tau/2)^2$. If we divide this expression by $\varepsilon$, will get the answer in the form $\tau = \varepsilon f(E,v_0)$ which is linear in $\varepsilon$. Commented Apr 2, 2014 at 16:27
• @Qmechanic, can you please recommend me a book where I find further discussion on this topic? Commented Jul 18, 2014 at 23:55
• @dapias: I assume you mean 1) the topic of regularization and 2) the topic of model building, i.e. what simplifying assumptions makes a good model. It is hard to recommend a few good concise books. Commented Jul 19, 2014 at 10:33

I think I got how to deal with the problem in a straightforward way, without the passage to a limit.

Let the phase space of the initial problem be the half plane

$$\{ \,(q,p) \; | \; q > 0\},$$

the wall is at $q=0$. In this phase space when the particle reaches the point $(0, p)$ it instantaneously teleports to the point $(0, -p)$. Particle trajectories are thus discontinuous.

The trick now is to glue the half plane into a cone, so that particle trajectories will become continuous. At the same time the particles are considered free on the entire trajectory, yielding the Hamiltonian being just a kinetic energy. Sorry for the lack of appropriate drawings, but I hope it is not so hard to imagine. The operation can be formally achieved with a non-canonical change of coordinates to $(r, \varphi)$:

$$\begin{cases} p = 2 r \sin \dfrac{\varphi}{2}, \\[.5em] q = 2 r \cos \dfrac{\varphi}{2}. \end{cases}$$

These are actually polar coordinates in a plane perpendicular to the cone axis. The symplectic form $dp \land dq$ transforms to $2r \, d \varphi \land dr$, or, in other words, the Poisson matrix becomes

$$\begin{bmatrix} 0 & \dfrac{1}{2r} \\ -\dfrac{1}{2r} & 0 \end{bmatrix}$$

The Hamiltonian (kinetic energy) is given by

$$H = \frac{2}{m} r^2 \sin^2 \frac{\varphi}{2}.$$

All of the above leads to the Hamiltonian flow of

$$X_H = \frac{1}{2m} r \sin \varphi \frac{\partial}{\partial r} - \frac{1 - \cos \varphi}{m} \frac{\partial}{\partial \varphi}$$

and the equations of motion of

$$\begin{cases} \dfrac{d r}{dt} = \dfrac{1}{2m} r \sin \varphi, \\[.5em] \dfrac{d \varphi}{dt} = -\dfrac{1}{m} (1 - \cos \varphi). \end{cases}$$

These are amenable to a relatively simple integration, free of the subtleties of special functions. As result one gets as a solution

$$\begin{gather} \cot \dfrac{\varphi(t)}{2} = C_1+ \dfrac{t}{m}, \\[.5em] r(t) = C_2 \sqrt{1 + \cot^2 \dfrac{\varphi(t)}{2}}, \end{gather}$$

which also can be obtained directly from the known solution in $(q,p)$ half-plane and coordinate transformation rules.

To sum up, the effect of the wall is accounted for via the change of the phase space topology. Thus, particles are considered free, with a Hamiltonian being a kinetic energy which is preserved. As far as I reckon the phase space is not a cotangent bundle of a configuration space anymore. If this is true along with all the above derivation, this represents probably the most simple case of a phase space that is not a cotangent bundle.

# Further investigation

Although at first I was guided by geometrical reasoning about the phase space, now I have been thinking about the coordinate transformation itself. I came up with another transform, which is much closer to the original coordinates:

$$\begin{cases} p = \operatorname{sgn} \! \left(\, \widetilde q \,\right) \, \widetilde p, \\[.5em] q = \left|\, \widetilde q \, \right| . \end{cases}$$

Here it is assumed that $\frac{d}{dx} \operatorname{sgn} x = 0$. Then $dp \land dq = d \, \widetilde p \land d \, \widetilde q$ and the Hamiltonian is the one of a free particle:

$$H = \frac{\; \widetilde p^{\, 2} \!}{2m}.$$

Like with trigonometric coordinate transformation there seems to be a need to choose the right branch for the transform $q \mapsto \widetilde q$, but since $q \geqslant 0$ there is no ambiguity.

• Distributions were invented for such problems, why go through all this trouble to avoid them? Commented Jun 8, 2014 at 20:59
• @auxsvr the problem is I do not see how they really solve the problem --- whenever I try to rigorously conduct the derivation I got stuck. And it seems to me that others hold the same opinion. I have commented on your answer on certain things I do not understand. Once again what is being sought --- Hamiltonian equations which can be solved as is, without resorting to any supplementary physical arguments. You see, Hamiltonian equations should be self contained, fully describing the system evolution, right? Commented Jun 9, 2014 at 5:47
• $\frac{d}{dx} \text{sgn}(x) = 2 \delta(x)$. Commented Jun 15, 2014 at 11:34
• I think auxsvr has a valid point, but your method is beautiful and highly ingenious! I loved reading this! Commented Jun 8, 2017 at 14:25

The problem here is not about a correct "physical" formalism. The problem is that we are seeking singular solutions. Let me explain (however, using Newtonian formalism, just because I'm more comfortable with it).

When talking about an initial value problem of an ODE of the form: $$\ddot{x}(t) = \frac{d}{dx}E(x(t))$$ one should have a uniform Lipschitz continuity of the right hand side to have a unique local classical solution according to the Picard-Lindelöf theorem. Clearly, in your case of a delta potential (i.e. hardwall) this is not the case since we have $$\ddot{x}(t)=-\delta^{'} (x(t))$$.

A common practice among mathematicians in such a case is to seek weak solutions (e.g. in the sense of Schwartz-Distribution theory). Instead of looking for a smooth solution $$x(t)$$, one looks for less regular solutions. This is in perfect accordance with the physical insight since we are expecting particles that have instantaneous change of momentum. To write the weak formulation of the problem we define $$\mathcal{D}(\mathbb{R})$$ to be the space of smooth functions with compact support and $$\mathcal{D}^{'}(\mathbb{R})$$ to be its dual space. The weak formulation reads: find $$x \in \mathcal{D}^{'}(\mathbb{R})$$ such that: $$\int \ddot{x} (s) \phi(s) \; ds = - \int \delta^{'} (x(s)) \phi(s) \; ds$$ where $$\phi \in \mathcal{D}(\mathbb{R})$$. Now, using integration by parts, one can relax the conditions on $$x(t)$$. A solution to the last equation is said to be a weak solution of our initial problem. However, here we have a problem resulting from the linear nature of Schwartz' distribution theory. The composition of a distribution $$\delta ^{'}$$ with an element $$x(t)$$ is only defined when $$x(t)$$ is smooth. In other words, we can not seek a weak solution $$x(t)$$ and at the same time consider the potential to be a distribution.

A proper and rigorous way to deal with this problem is using a nonlinear theory of distributions, where the composition of a distribution is defined, even for less regular objects. This is possible in Colombeau's algebra of tempered generalised functions.

In "Geometric Theory of Generalised Functions with Applications to General Relativity" by Michael Grosser et al. Newton's equation with a delta potential was proven to have a unique solution in Colombeau's algebra. This solution is parametrised by a regularisation parameter $$\varepsilon$$. Moreover, it's been proven that when taking the limit of the solution as $$\varepsilon \to 0$$ it converges to the solution one expects in physics (i.e. the trajectory experiences a reflection at the potential).

Note that this is very similar to the regularisation procedure described in one of the answers. However, in this algebra the regularisation is made more rigorous and consistent with distribution theory.

I will yet add another formalism:

Lets start with the hamiltonian form of Hamilton's Principle. Let $c: \mathbb R \longrightarrow T^*Q; t\mapsto (q(t),p(t))$ be the trajectory of a particle in the phase space of the configuration space $Q$, we define a subset of $Q$, $C$, where no contac occur between the particles, and $\partial C$ is te set of ponts wer contac has occure, but not penetration. We assume that $c$ will intersect $\partial C$ at time $t_c$, that is $q(t_c)\in \partial C$. $c$ is smooth out of $t_c$. Lets now derive the equations of motin with variations. We start with the action in hamiltoninan formalism with Lagrangian $L = \dot q p -H(q,p)$:

$$S[c] = \int_0^T L(q,\dot q, p)dt = \int_0^{t^\lambda} L(q,\dot q, p)dt + \int_{t^\lambda} ^T L(q,\dot q, p)dt ,$$ and with the variation on the path denoted by $c^\lambda$ and the variation on the collision instant by $t^\lambda$, if we set $\delta S =0$, by the Leibniz integration rule: $$\frac{d}{d\lambda}S[c^\lambda] = \int^{t^\lambda}_0 \left[\frac{\partial L}{\partial c}\delta c +\frac{\partial L}{\partial \dot c}\delta \dot c \right] dt + \int^{T}_{t^\lambda} \left[\frac{\partial L}{\partial c}\delta c +\frac{\partial L}{\partial \dot c}\delta \dot c \right] dt + \left.L\delta t^\lambda\right|_{t^\lambda_-} -L\delta$$

t^\lambda\right|{t^\lambda+} $$We integrate by parts as usual and require that \delta q(0)=\delta q(T)=0, there is no need to \delta p to vanish as it is not integrated by parts, we get:$$ \int^{t^\lambda}_0 \left[\frac{\partial L}{\partial c}\delta c -\frac{d}{dt}\frac{\partial L}{\partial \dot c}\delta c \right] dt + \int^{T}_{^\lambda} \left[\frac{\partial L}{\partial c}\delta c -\frac{d}{dt}\frac{\partial L}{\partial \dot c}\delta c \right] dt - \left[\frac{\partial L}{\partial \dot c}\delta c + L\delta t^\lambda\right]_{t_-^\lambda}^{t^\lambda_+}, $$as we require the variations to vanish at all time, we can focus only in the third term:$$ \left[\frac{\partial L}{\partial \dot c}\delta c + L\delta t_c\right]_{t_c^-}^{t_c^+} = \left[p \delta q + (\dot q p - H)\delta t_c\right]_{t_c^-}^{t_c^+}. $$We recall that q(t_c)\in \partial C, so the variation \delta q(t_c) + \dot q(t_c)\delta t_c \in T_{q(t_c)}Q and we get:$$ \left[p \delta q + (\dot q p - H)\delta t^\lambda\right]_{t_-^\lambda}^{t_+^\lambda} = \left[p(\delta q + \dot q p\delta t^\lambda) - H \delta t^\lambda \right]_{t_-^\lambda}^{t_+^\lambda} =0, $$considering the variations arbitrary, we get conservation of energy, H has to be constant and we get that the change of momentum p(t_+) -p(t_-) has to be normal to T_{q(t_c)}Q, and the conservation of energy fixes the magnitude. The major cave beat is that it is necessary to interpret which direction the particle exits the collision, because it is not possible to extract the information from the equations Sources: http://thesis.library.caltech.edu/1934/1/01thesis.pdf http://people.math.sfu.ca/~van/papers/40603.pdf Try using U(x) = \Theta(x), \Theta(x) \equiv \begin{cases} 1, & x>0\\ 0, & x<0\end{cases} . The force becomes \dot{p} = - \delta(x), as anticipated. Because the force is infinite with direction opposite that of the momentum, the particle cannot cross the plane x=0. Another way to see that the particle won't cross x=0 is to integrate the equation of motion over x just before and just after x=0:$$ -1 = m\int_{0_-}^{0_+} \dot{v} dx = \frac{m}{2} (v^2(0_+)- v^2(0_-)), from which we obtain $v(0_+) = 0$, because the energy is $1$.

• Could you please elaborate on how you have integrated $\int_{0_-}^{0^+} \dot \nu \, dx$ ? This is not so obvious. Commented Jun 9, 2014 at 4:51
• Oh, I got it, but still it would be nicer to be written. Moving along the derivation, as far as I understand you show that momentum when $x > 0$ is $0$ but it does not justify that is should reverse. Andy why is the energy is $1$ as you said? It should be $p^2/(2m)$ in any case. Moreover, with the potential fixed being equal $1$ does not prohibit the particle to enter the region if its kinetic energy is greater then $1$. Commented Jun 9, 2014 at 5:40
• I used a finite potential, because it is the more general case, provided that $U(x) = U_{\text{max}} \Theta(x) \geq p^2/2m$, which can be considered to be the definition of the problem, otherwise the particle will cross the potential. Thus, the last equation becomes $0 = \frac{m}{2} v^2(0_+) + [ U_{\text{max}} - \frac{m}{2} v^2(0_-)]$, which implies that $v(0_+) =0$ and $U_{\text{max}} = \frac{m}{2} v^2(0_-)$. Now you have a particle with nonzero momentum only before the barrier and a force pointing in the negative direction. This is enough to conclude that it will reverse motion. Commented Jun 9, 2014 at 7:55
• If you require only $U_\text{max}$ being greater or equal than the initial kinetic energy, then $0=\frac{m}{2} \nu^2(0_+) + [U_\text{max} − \frac{m}{2} \nu^2(0_+)]$ may not have solutions at all with $\nu(0_+) \in \mathbb R$. It will have a solution $\nu(0_+) = 0$ only if $U_\text{max}$ equals to a kinetic energy of an incident particle, but that makes it dependent on the particle momentum. Commented Jun 9, 2014 at 8:33
• Yes, but this is part of the definition of the problem, $U_\text{max} \geq p^2/2m$, otherwise the particle will cross the potential. The previous argument establishes that $\exists U_\text{max}$ such that the particle does not cross the barrier. Since for $U_\text{max} >p^2/2m$ the equation may not have a real solution, we resort to checking whether the particle crosses the barrier for increasing values of $U_\text{max}$ with $U_\text{max} < p^2/2m$, which implies that $v(0_+)$ decreases as $U_\text{max}$ increases, and, since higher potential means greater opposing force, Commented Jun 9, 2014 at 9:34