Denote the length distribution of one’s hair to be f(l,t)f(l,t). This means that, at time tt, the number of hairs within the length range from ll to l+dll+\mathrm dl is Nf ⁣(l,t)dlNf\!\left(l,t\right)\mathrm dl, where NN is the total number of hairs.

Each hair grows at constant speed vv. However, they cannot grow indefinitely because there is a probability of λ\lambda per unit time for a hair to be lost naturally (this is the same assumption as the exponential decay). After a hair is lost, it restarts growing from zero length.

Suppose that at t=0t=0 you have got a haircut so that the hair length distribution becomes

f(l,0)=f0(l).f(l,0)=f_0(l). (1)(1)

Then, how does the distribution evolve with time?

Because ff is a distribution function, There is a normalization restriction on ff:

∫0∞f ⁣(l,t)dl=1,t≥0.\int_0^\infty f\!\left(l,t\right)\mathrm dl=1,\quad t\ge0. (2)(2)

This normalization condition also applies to the initial condition (t=0t=0). This means that the function f0f_0 also satisfies the normalization restriction (Equation 2)

∫0∞f0 ⁣(l)dl=1.\int_0^\infty f_0\!\left(l\right)\mathrm dl=1. (3)(3)

Because of the natural loss of hair, only a 1−λ dt1-\lambda\,\mathrm dt portion of hair will survive dt\mathrm dt. According to this, we can construct the following equation:

f ⁣(l+v dt,t+dt)=(1−λ dt)f ⁣(l,t).f\!\left(l+v\,\mathrm dt,t+\mathrm dt\right)= \left(1-\lambda\,\mathrm dt\right)f\!\left(l,t\right).

This equation can be reduced to a first-order linear PDE:

v∂f∂l+∂f∂t+λf=0.v\frac{\partial f}{\partial l}+\frac{\partial f}{\partial t}+\lambda f=0. (4)(4)

Define θ≔l−vt\theta\coloneqq l-vt, and define a new function

g ⁣(θ,t)≔f ⁣(θ+vt,t).g\!\left(\theta,t\right)\coloneqq f\!\left(\theta+vt,t\right).

Then,

∂g∂t=v∂f∂l+∂f∂t.\frac{\partial g}{\partial t}=v\frac{\partial f}{\partial l}+\frac{\partial f}{\partial t}.

Then, Equation 4 can be reduced to

∂g∂t+λg=0.\frac{\partial g}{\partial t}+\lambda g=0.

The solution is

g ⁣(θ,t)=Φ ⁣(θ)e−λt,g\!\left(\theta,t\right)=\Phi\!\left(\theta\right)\mathrm e^{-\lambda t},

where Φ\Phi is an arbitrary function defined on R\mathbb R. Therefore, the general solution to Equation 4 is

f ⁣(l,t)=Φ ⁣(l−vt)e−λt.f\!\left(l,t\right)=\Phi\!\left(l-vt\right)\mathrm e^{-\lambda t}. (5)(5)

By utilizing Equation 1, we can find Φ(θ)\Phi(\theta) for θ>0\theta>0. Substitute t=0t=0 into Equation 5 and compare with Equation 1, and we have

Φ(θ>0)=f0(θ).\Phi(\theta>0)=f_0(\theta). (6)(6)

This only gives Φ(θ)\Phi(\theta) for θ>0\theta>0 because f0f_0 is not defined on negative numbers. The rest of Φ\Phi, however, may be deduced from Equation 2.

Substitute Equation 5 into Equation 2, and we have

1=∫0∞Φ ⁣(l−vt)e−λt dl=e−λt(∫0vtΦ ⁣(l−vt)dl+∫vt∞Φ ⁣(l−vt)dl)=e−λt(∫−vt0Φ ⁣(θ)dθ+∫0∞f0 ⁣(θ)dθ)=e−λt(1−∫0−vtΦ ⁣(θ)dθ).\begin{align*} 1&=\int_0^\infty\Phi\!\left(l-vt\right)\mathrm e^{-\lambda t}\,\mathrm dl\\ &=\mathrm e^{-\lambda t}\left( \int_0^{vt}\Phi\!\left(l-vt\right)\mathrm dl +\int_{vt}^\infty\Phi\!\left(l-vt\right)\mathrm dl \right)\\ &=\mathrm e^{-\lambda t}\left( \int_{-vt}^0\Phi\!\left(\theta\right)\mathrm d\theta +\int_0^\infty f_0\!\left(\theta\right)\mathrm d\theta \right)\\ &=\mathrm e^{-\lambda t}\left(1-\int_0^{-vt}\Phi\!\left(\theta\right)\mathrm d\theta\right). \end{align*}

Therefore, we have the integral of Φ\Phi on negative intervals:

∫0−vtΦ ⁣(θ)dθ=1−eλt.\int_0^{-vt}\Phi\!\left(\theta\right)\mathrm d\theta =1-\mathrm e^{\lambda t}.

Find the derivative of both sides of the equation w.r.t. tt, and we have

−vΦ ⁣(−vt)=−λeλt.-v\Phi\!\left(-vt\right)=-\lambda\mathrm e^{\lambda t}.

In other words,

Φ ⁣(θ<0)=λve−λvθ.\Phi\!\left(\theta<0\right)=\frac\lambda v\mathrm e^{-\frac\lambda v\theta}. (7)(7)

Combining Equation 6 and 7 and substituting back to Equation 5, we can find the special solution to Equation 4 subject to restrictions Equation 2 and 1:

f(l,t)={λve−λvl,0<l<vt,e−λtf0(l−vt),l>vt.f(l,t)=\begin{cases} \frac\lambda v\mathrm e^{-\frac\lambda vl},&0<l<vt,\\ \mathrm e^{-\lambda t}f_0(l-vt),&l>vt. \end{cases} (8)(8)

This is our final answer.

This result is interesting in that any distribution will finally evolve into an exponential distribution with the rate parameter being λv\frac\lambda v as t→∞t\to\infty no matter what the initial distribution is:

f(l,∞)=λve−λvl.f(l,\infty)=\frac\lambda v\mathrm e^{-\frac\lambda vl}.

This distribution is the stationary solution to Equation 4. This is actually a normal behavior for first-order PDEs. For example, the thermal equilibrium state is the stationary solution to the heat equation, and any other solution approaches to the stationary solution over time.

This behavior can explain why human body hair tends to grow to only a certain length instead of being indefinitely long. You may try shaving your leg hair and wait for some weeks. You can observe that they grow to approximately the original length but not any longer. It is similar for your hair (on top of your head), but λ\lambda of hair is so small that it can hardly reach its terminal length if you get haircuts regularly.

Another thing to note is that this may explain a phenomenon that we may observe: the longer your hair is, the more slowly it grows, and your hair no longer seems to grow when it reaches a certain length. If the length of hair that we observe is actually the mean length of the hair, then it is

μ(t)=∫0∞lf ⁣(l,t)dl=(μ0−vλ)e−λt+vλ,\mu(t)=\int_0^\infty lf\!\left(l,t\right)\mathrm dl =\left(\mu_0-\frac v\lambda\right)\mathrm e^{-\lambda t}+\frac v\lambda,

where μ0\mu_0 is the mean of the distribution f0f_0. It can be seen that the growth rate of the mean length of hair varies exponentially.