Best I could determine after searching for a couple of hours, this function was first introduced in the following publication, a scan of which is made available by the Bibliothèque nationale de France:
P. Langevin, "Magnétisme et théorie des électrons," Annales de chimie et de physique, 8:5 (1905), 70-127.
On page 117 one finds:
Le moment magnétique total des $\mathrm{N}$ molécules est évidemment
dirigé parallèlement au champ et est égal à la somme des projections
sur cette direction des moments composants. Pour l'unité de volume
supposée contenir $\mathrm{N}$ molécules, ce moment résultant
représente l'intensité d'aimantation $\mathrm{I}$ :
$$\mathrm{I}=\int{\mathrm{M}\cos \alpha}{\ dn} = \int_{-1}^{+1}{2\pi\mathrm{MK\ }x\ e^{ax}\ {dx}}.$$
Or
$$\int_{-1}^{+1}{x\ }{e^{ax}}{\ dx} = 2 \left(\frac{\mathrm{Ch\ } {a}}{a}-\frac{\mathrm{Sh\ }a}{a^{2}}\right),$$
d'où $$\mathrm{I}=\mathrm{MN}\left(\frac{\mathrm{Ch\ }a}{\mathrm{Sh\
}{a}}-\frac{1}{a}\right).$$
Here $\mathrm{Ch}$ represents the hyperbolic cosine and $\mathrm{Sh}$ the hyperbolic sine, thus $\mathrm{I}=\mathrm {MN}\left(\mathrm{coth}(a) - \frac{1}{a}\right) = \mathrm{MN}\mathcal{L}(a).$