# Structure of the probability mass function of the Poisson distribution of
order $k$

S. R. Mane<EMAIL_ADDRESS>Convergent Computing Inc., P. O. Box 561,
Shoreham, NY 11786, USA

###### Abstract

The Poisson distribution of order $k$ is a special case of a compound Poisson
distribution. For $k=1$ it is the standard Poisson distribution. Although its
probability mass function (pmf) is known, what is lacking is a visual
interpretation, which a sum over terms with factorial denominators does not
supply. Unlike the standard Poisson distribution, the Poisson distribution of
order $k$ can display a maximum of four peaks simultaneously, as a function of
two parameters: the order $k$ and the rate parameter $\lambda$. This note
characterizes the shape of the pmf of the Poisson distribution of order $k$.
The pmf can be partitioned into a single point at $n=0$, an increasing
sequence for $n\in[1,k]$ and a mountain range for $n>k$ (explained in the
text). The “parameter space” of the pmf is mapped out and the significance of
each domain is explained, in particular the change in behavior of the pmf as a
domain boundary is crossed. A simple analogy (admittedly unrelated) is that of
the discriminant of a quadratic with real coefficients: its domains
characterize the nature of the roots (real or complex), and the domain
boundary signifies the presence of a repeated root. Something similar happens
with the pmf of the Poisson distribution of order $k$. As an application, this
note explains the mode structure of the Poisson distribution of order $k$.
Improvements to various inequalities are also derived (sharper bounds, etc.).
New conjectured upper and lower bounds for the median and the mode are also
proposed.

###### keywords:

Poisson distribution of order $k$ , probability mass function , median , mode
, Compound Poisson distribution , discrete distribution

###### MSC:

[2020] 60E05 , 39B05 , 11B37 , 05-08

††journal: (internal report CC23-6)

##  1  Introduction

In two recent notes [1, 2], the author presented numerical results for the
Poisson distribution of order $k$ [3]. It is a variant (or extension) of the
well-known Poisson distribution. We begin with its formal definition.

######  Definition 1.1.

The Poisson distribution of order $k$ (where $k\geq 1$ is an integer) and
parameter $\lambda>0$ is an integer-valued statistical distribution with the
probability mass function (pmf)

$f_{k}(n;\lambda)=e^{-k\lambda}\sum_{n_{1}+2n_{2}+\dots+kn_{k}=n}\frac{\lambda^{n_{1}+\dots+n_{k}}}{n_{1}!\dots
n_{k}!}\,,\qquad n=0,1,2\dots$ (1.1)

For $k=1$ it is the standard Poisson distribution. The Poisson distribution of
order $k$ is a special case of the compound Poisson distribution introduced by
Adelson [4]. Although exact expressions for the mean and variance of the
Poisson distribution of order $k$ are known [5], exact results for its median
and mode are difficult to obtain.

What is lacking is a visual interpretation of the pmf, which a formal
mathematical sum such as eq. (1.1) does not supply. Unlike the standard
Poisson distribution, the Poisson distribution of order $k$ can display a
maximum of four peaks simultaneously, as a function of two parameters: the
order $k$ and the rate parameter $\lambda$. This note characterizes the shape
of the pmf of the Poisson distribution of order $k$. The “parameter space” of
the pmf is mapped out and the significance of each domain is explained, in
particular the change in behavior of the pmf as a domain boundary is crossed.
A simple analogy (admittedly unrelated) is that of the discriminant of a
quadratic with real coefficients: its domains characterize the nature of the
roots (real or complex), and the domain boundary signifies the presence of a
repeated root. Something similar happens with the pmf of the Poisson
distribution of order $k$. As an application, this note explains the mode
structure of the Poisson distribution of order $k$. Improvements to various
inequalities are derived (sharper bounds, etc.). New conjectured upper and
lower bounds for the median and the mode are also proposed.

The structure of this paper is as follows. Sec. 2 presents basic definitions
and notation employed in this note. Improvements to some published
inequalities are also presented (sharper bounds, etc.) Sec. 3 quantifies the
structure of the pmf of the Poisson distribution of order $k$. Secs. 4 and 5
present new conjectures for upper and lower bounds for the median and mode,
respectively. Sec. 6 concludes.

##  2 Basic notation and definitions

For later reference we define the parameter $\kappa=k(k+1)/2$. We denote the
mean by $\mu$, the median by $\nu$ and the mode by $m$ (with pertinent
subscripts, etc. to denote the dependence on $k$ and $\lambda$, see below).
Philippou [5] derived that the mean is $\mu_{k}(\lambda)=\kappa\lambda$ and
the variance is $\sigma_{k}^{2}(\lambda)=\frac{1}{6}k(k+1)(2k+1)\lambda$. For
the median, we follow the exposition in [6]: if $Y_{k,\lambda}$ is a random
variable which is Poisson distributed with order $k$ and parameter $\lambda$,
the median is defined as the smallest integer $\nu$ such that
$P(Y_{k,\lambda}\leq\nu)\geq\frac{1}{2}$. With this definition, the median is
unique and is always an integer. The mode is defined as the location(s) of the
global maximum of the probability mass function. It is known that the mode may
not be unique. For the standard Poisson distribution with parameter $\lambda$,
the mode equals $\lfloor\lambda\rfloor$ if $\lambda\not\in\mathbb{N}$, but
both $\lambda-1$ and $\lambda$ are modes if $\lambda\in\mathbb{N}$. We adopt
the following notation and definitions from [7].

  1. 1.

We work with $h_{k}(n;\lambda)=e^{k\lambda}f_{k}(n;\lambda)$ (see [7]) and
refer to it as the “scaled pmf” below. Observe from eq. (1.1) that
$h_{k}(n;\lambda)$ is a polynomial in $\lambda$ with all positive
coefficients. It has degree $n$ and for $n>0$ it has no constant term (also
$h_{k}(0;\lambda)=1$ and $h_{k}(1;\lambda)=\lambda$ for all $k\geq 1$). Hence
for fixed $k\geq 1$ and $n>0$, $h_{k}(n;0)=0$ and $h_{k}(n;\lambda)$ is a
strictly increasing function of $\lambda$ for $\lambda>0$.

  2. 2.

The parameter $r_{k}$ is defined as the positive root of the equation
$h_{k}(k;\lambda)=1$. It was shown in [7] that $r_{k}$ is unique and
$0<r_{k}<1$.

  3. 3.

It was proved (Lemma 1 in [7]) that for fixed $k\geq 2$ and $\lambda>0$, the
sequence $\\{h_{k}(n;\lambda),n=1,\dots,k\\}$ is strictly increasing, i.e.
$h_{k}(n-1;\lambda)<h_{k}(n;\lambda)$ for $n=2,\dots,k$. Hence only the last
index $k$ can be a mode. An integer in the interval $[1,k-1]$ can never be a
mode of the Poisson distribution of order $k$.

  4. 4.

It was proved (Lemma 3 in [7]) that for fixed $k\geq 2$ and $0<\lambda\leq
r_{k}$, then $h_{k}(k;\lambda)>h_{k}(k+1;\lambda)$. This makes
$h_{k}(k;\lambda)$ a local maximum in the histogram of the pmf, for
sufficiently small values of $\lambda$. We shall see this below, when plotting
graphs of the histogram of the pmf. Note that the condition $0<\lambda\leq
r_{k}$ is sufficient but not necessary to attain
$h_{k}(k;\lambda)>h_{k}(k+1;\lambda)$.

By the term “double mode” we mean the distribution is bimodal, with joint
modes at $m_{1}$ and $m_{2}$. For $k=1$, the standard Poisson distribution,
the integers $m_{1}$ and $m_{2}$ are always consecutive integers, but for
$k\geq 2$ this need not be so. Kwon and Philippou [7] tabulated a list of
double modes for $k=2,3,4$ and $0<\lambda\leq 2$. It was shown in [1] that for
any $k\geq 2$, the Poisson distribution of order $k$ has a denumerable
infinity of double modes, consisting of pairs of consecutive integers. A major
topic of this note is to characterize the mode structure for $k\geq 2$,
including double modes with non-consecutive integers. The existence of three
or more joint modes is an open question. This note presents numerical evidence
that the Poisson distribution of order $k$ does not have three or more joint
modes.

The term “first double mode” signifies the first time (smallest value of
$\lambda$) that the Poisson distribution of order $k$ has a double mode. The
mode values in this case are $0$ and $m>0$. The following notation was
introduced in [2] for the first double mode. The nonzero mode value was
denoted by $\hat{m}_{k}$ and the corresponding value of $\lambda$ was denoted
by $\hat{\lambda}_{k}$. In terms of this notation, Kwon and Philippou [7]
showed that $\hat{m}_{k}=k$ for $k=2,3,4$. It was shown in [1] that
$\hat{m}_{k}=k$ for $k=2,\dots,14$ and $\hat{m}_{k}>k$ for $k\geq 15$. It was
proved in [2] that $\hat{m}_{k}\geq k$ for all $k\geq 2$, but an exact formula
for $k\geq 15$ is not known to date. Nor is it proved that $\hat{m}_{k}$
increases strictly with $k$ (for $k\geq 15$) although numerical tests up to at
least $k=10^{4}$ have not found any exceptions [2].

We take the opportunity here to present improvements to various published
inequalities (sharper bounds, etc.). The following inequality was derived in
[2]

$1/\kappa\leq\hat{\lambda}_{k}\leq r_{k}<1\,.$ (2.1)

Recall that it was proved in [7] that for fixed $k\geq 2$ and $\lambda>0$, the
sequence $\\{h_{k}(n;\lambda),n=1,\dots,k\\}$ is strictly increasing. It
follows that the location of the first double mode is not less than $k$, i.e.
$\hat{m}_{k}\geq k$. Next, the mode is bounded by the floor of the mean
(Theorem 2.1 in [8], recall eq. (5.2)) and Philippou [5] showed that the value
of the mean is $\kappa\lambda$. Hence

$\hat{m}_{k}\leq\kappa\hat{\lambda}_{k}\,.$ (2.2)

We can deduce two inequalities from this information.

######  Proposition 2.1.

Using eq. (2.2) and $\hat{\lambda}_{k}<1$, it follows that for all $k\geq 2$,

$k\leq\hat{m}_{k}<\kappa\,.$ (2.3)

######  Proposition 2.2.

Using eq. (2.2) and $\hat{m}_{k}\geq k$, we deduce
$\kappa\hat{\lambda}_{k}\geq k$. Solving for $\hat{\lambda}_{k}$ yields the
inequality

$\hat{\lambda}_{k}\geq\frac{k}{\kappa}=\frac{2}{k+1}\,.$ (2.4)

######  Remark 2.3.

Using eq. (2.4), we improve the inequalities in eq. (2.1) as follows

$\frac{2}{k+1}\leq\hat{\lambda}_{k}\leq r_{k}<1\,.$ (2.5)

######  Remark 2.4.

Philippou [9] showed that the Poisson distribution of order $k$ has a unique
mode of zero if $\lambda<1/\kappa=2/(k(k+1))$. Using eq. (2.4), we improve
this bound to say the mode is uniquely zero if

$\lambda<\frac{2}{k+1}\,.$ (2.6)

This is a sufficient but not necessary condition.

######  Remark 2.5.

We can now prove the conjecture in [1] that if the median is zero then the
mode is also zero.

###### Proof.

The median is zero if and only if $\lambda\leq(\ln 2)/k$ (proved in [1]).
Observe that for all $k\geq 1$,

$\frac{\ln 2}{k}<\frac{2}{k+1}\,.$ (2.7)

Hence if $\lambda\leq(\ln 2)/k$ (median is zero) then $\lambda<2/(k+1)$ and
from eq. (2.6) the mode is zero. ∎

##  3 Structure of the probability mass function

###  3.1 General remarks

We shall plot graphs to investigate the structure of the Poisson distribution
of order $k$. To fix ideas, we plot the value of the scaled pmf
$h_{k}(n;\lambda)$, as opposed to the true pmf $f_{k}(n;\lambda)$ in eq.
(1.1). The prefactor $e^{-k\lambda}$ does not alter the shape of the histogram
and is therefore irrelevant to the analysis below. To demonstrate, the scaled
pmf of the standard Poisson distribution is plotted in Fig. 1.

  1. 1.

For $\lambda<1$, e.g. $\lambda=0.8$, the unique mode is zero and the pmf
decreases monotonically.

  2. 2.

When $\lambda=1$, the first double mode is attained, at $0$ and $1$.

  3. 3.

As the value of $\lambda$ increases further, the pmf displays a single peak.
Its location shifts to the right as the value of $\lambda$ increases. As the
value of $\lambda$ increases, the location of the mode increases in unit
steps.

  4. 4.

If the value of $\lambda$ is an integer, e.g. $\lambda=4$, the peak has a flat
top and the distribution is bimodal, with joint modes at $\lambda-1$ and
$\lambda$.

  5. 5.

If the value of $\lambda$ is not an integer, e.g. $\lambda=4.2$, the mode is
unique, with value $\lfloor\lambda\rfloor$.

By contrast, Fig. 2 displays a plot of the scaled pmf of the Poisson
distribution of order $50$, for $\lambda\simeq 0.10194$. There are four peaks
in the histogram, at (i) $n=0$ (height$\;=1$), (ii) $n=50$ (height$\;\simeq
0.6698$), which is a local maximum and not a mode, (iii) $n=98$
(height$\;\simeq 0.98358$), which is also a local maximum and not a mode, (iv)
$n=113$ (height$\;=1$), which is a double mode along with $0$. As the value of
$\lambda$ changes, the relative heights of the peaks change, and the mode
structure will vary. For example, we know that for sufficiently small values
$0<\lambda<2/(k+1)$, the unique mode is $0$. Our goal in this section is to
understand the evolution of the structure of the scaled pmf of the Poisson
distribution of order $k$.

As complicated as it appears, Fig. 2 does contain useful clues as to how to
parameterize the scaled pmf. First, the point at $n=0$ always has a height of
$1$ and is an invariant: it is the same for all $k$ and $\lambda$. Next, it
was mentioned earlier that the points for $1\leq n\leq k$ form a strictly
increasing sequence (this was proved in Lemma 1 in [7]) and that sequence is
visible in Fig. 2. It was also mentioned earlier that for sufficiently small
$\lambda$, the point at $n=k$ is a local maximum (proved in Lemma 3 in [7])
and this fact is also visible in Fig. 2. The remaining points $n>k$ exhibit a
more or less “mountain range” appearance, similar to the standard Poisson
distribution, except there are two peaks. We shall refer to the region $n>k$
as the “mountain range” region and speak of the “left peak” and the “right
peak” below. However, be advised that the mountain range does not always have
two peaks: that is why we need a “parameter map” to classify matters.

We wish to compose a “parameter map” to characterize the behavior of the
scaled pmf of the Poisson distribution of order $k$ as a function of the order
$k$ and the rate parameter $\lambda$. For the dependence on $k$, there are
four cases: (i) $k\in[2,3]$, (ii) $k\in[4,14]$, (iii) $k\in[15,41]$, (iv)
$k\geq 42$. In each case, we increase the value of $\lambda$ continuously from
zero and examine the structure of the scaled pmf, and the resulting
consequences.

###  3.2 Case $k\in[2,3]$

We select $k=3$ as a representative example, because for $k=2$ the “increasing
sequence” $n=1,\dots,k$ contains only two points and is not informative. The
scaled pmf of the Poisson distribution of order $3$ is plotted in Fig. 3 for
selected values of $\lambda$, increasing from top to bottom.

  1. 1.

In the top panel, $\lambda=0.4$. This value is so small that the mountain
range region ($n>k$) is monotonically decreasing. The increasing sequence for
$n=1,\dots,k$ is visible, as is the local maximum at $n=k$, but it is not a
global maximum. The unique mode is zero.

  2. 2.

As the value of $\lambda$ increases, the height of the local maximum at $n=k$
reaches $1$ and a double mode is attained. The joint modes are at $0$ and $k$.
This is displayed in the second panel, where $\lambda\simeq 0.601679$. The
mountain range region ($n>k$) is monotonically decreasing, although one can
discern that a peak is forming.

  3. 3.

When the value of $\lambda$ increases, the unique mode jumps from $0$ to $k$,
i.e. an increase of more than one unit. The “first double mode” represents a
boundary between two domains in the parameter space.

  4. 4.

As the value of $\lambda$ increases further, the mountain range develops a
single peak, and its height rises to equal that of the point at $n=k$. This is
displayed in the third panel, where $\lambda\simeq 0.9962$. The joint modes
are at $k$ and $k+2$ (for both cases $k=2$ and $k=3$).

  5. 5.

When the value of $\lambda$ increases, the unique mode jumps from $k$ to
$k+2$, i.e. an increase of more than one unit. This is the “second double
mode” and is a boundary between two domains in the parameter space.

  6. 6.

As the value of $\lambda$ increases further, the height of the peak of the
mountain range exceeds that of the point at $n=k$. The location of the single
peak shifts rightwards and its height increases. The mode is determined solely
by the location of the peak of the mountain range and increases in unit steps
and takes all integer values $\geq k+2$. This is displayed in the fourth
panel, where $\lambda=1.02$. The scaled pmf has a unique mode given by the
height of the single mountain peak. The point at $n=k$ is still a local
maximum but plays no further role to determine the mode.

  7. 7.

For a discrete (denumerably infinite) set of values of $\lambda$, the mountain
peak has a flat top and the distribution is bimodal. The joint modes consist
of consecutive integers. This is displayed in the fifth panel, where
$\lambda\simeq 1.4293$. There is a double mode, at $7$ and $8$. It is almost a
triple mode, but the height of the point at $n=6$ is a little lower. The point
at $n=k$ is no longer a local maximum.

###  3.3 Case $k\in[4,14]$

We select $k=10$ as a representative example. The scaled pmf of the Poisson
distribution of order $10$ is plotted in Fig. 4 for selected values of
$\lambda$, increasing from top to bottom.

  1. 1.

In the top panel, $\lambda=0.2$. This value is so small that the mountain
range region ($n>k$) is monotonically decreasing. The increasing sequence for
$n=1,\dots,k$ is visible, as is the local maximum at $n=k$, but it is not a
global maximum. The unique mode is zero.

  2. 2.

As the value of $\lambda$ increases, the height of the local maximum at $n=k$
reaches $1$ and a double mode is attained. The joint modes are at $0$ and $k$.
This is displayed in the second panel, where $\lambda\simeq 0.31713$. The
mountain range region ($n>k$) exhibits a left peak, but its height is less
than $1$ and is a local but not global maximum.

  3. 3.

When the value of $\lambda$ increases, the unique mode jumps from $0$ to $k$,
i.e. an increase of more than one unit. The “first double mode” represents a
boundary between two domains in the parameter space.

  4. 4.

As the value of $\lambda$ increases further, the height of the left peak rises
and equals that of the point at $n=k$. This is displayed in the third panel,
where $\lambda\simeq 0.36189$. The joint modes are at $k$ and the location of
the left mountain peak, say $m_{\rm left}$. The value of $m_{\rm left}$
depends on $\lambda$ but is always at least $k+2$. No explicit formula is yet
known for $m_{\rm left}$ when a double mode is attained with the point at
$n=k$, although the values can be tabulated for all $k\in[4,14]$.

  5. 5.

When the value of $\lambda$ increases, the unique mode jumps from $k$ to
$m_{\rm left}$, i.e. an increase of more than one unit. This is the “second
double mode” and is a boundary between two domains in the parameter space.

  6. 6.

As the value of $\lambda$ increases further, the mode is determined by the
location of the left mountain peak. The value of the mode increases in unit
steps and there are double modes consisting of pairs of consecutive integers.
However, the mountain range develops a second (right) peak and its height
rises faster than that of the left peak and it catches up with the height of
the left peak.

  7. 7.

As the value of $\lambda$ increases further, a double mode is attained, where
the heights of the two mountain peaks are equal. This is displayed in the
fourth panel, where $\lambda\simeq 0.472694$. We denote the locations of the
two mountain peaks by $m_{\rm left}$ and $m_{\rm right}$, respectively. Very
significant: the value of $m_{\rm left}$ is larger than that in the previous
panel, i.e. the value of $m_{\rm left}$ increases with $\lambda$. The point at
$n=k$ is still a local maximum but plays no further role to determine the
mode. Observe also that there is a local minimum between the two mountain
peaks, i.e. they are separated by more than one unit and are not consecutive
integers.

  8. 8.

When the value of $\lambda$ increases, the unique mode jumps from $m_{\rm
left}$ to $m_{\rm right}$, which is an increase of more than one unit.
However, this is not the “third double mode” because the value of $m_{\rm
left}$ increased between the third and fourth panels, i.e. there were double
modes consisting of pairs of consecutive integers. We may refer to this as the
“third mode jump” (by which is meant an increase of the mode by more than one
unit) and is a boundary between two domains in the parameter space. With this
terminology, the first double mode is the first mode jump and the second
double mode is the second mode jump.

  9. 9.

As the value of $\lambda$ increases further, the height of the right mountain
peak exceeds that of the left. The location of the right peak shifts
rightwards and its height increases. The mode is determined solely by the
location of the right mountain peak and increases in unit steps and takes all
integer values from its value at the third mode jump upwards.

  10. 10.

For a discrete (denumerably infinite) set of values of $\lambda$, the right
mountain peak has a flat top and the distribution is bimodal. The joint modes
consist of consecutive integers. This is displayed in the fifth panel, where
$\lambda\simeq 0.5119$. There is a double mode at $24$ and $25$. The point at
$n=k$ is still a local maximum but is no longer relevant to determine the
mode.

###  3.4 Case $k\in[15,41]$

We select $k=20$ as a representative example. The scaled pmf of the Poisson
distribution of order $20$ is plotted in Fig. 5 for selected values of
$\lambda$, increasing from top to bottom.

  1. 1.

In the top panel, $\lambda=0.1$. This value is so small that the mountain
range region ($n>k$) is monotonically decreasing. The increasing sequence for
$n=1,\dots,k$ is visible, as is the local maximum at $n=k$, but it is not a
global maximum. The unique mode is zero.

  2. 2.

As the value of $\lambda$ increases, the mountain range develops a left peak
and its height rises and catches up with that of the point at $n=k$. This is
displayed in the second panel, where $\lambda\simeq 0.1899$. Note that both
heights are less than $1$, hence they are not modes.

  3. 3.

As the value of $\lambda$ increases further, the height of the left mountain
peak rises faster than that of the point at $n=k$ and it equals $1$ before the
point at $n=k$ does so. This is displayed in the third panel, where
$\lambda\simeq 0.20333$. Hence the first double mode is not at $0$ and $k$.
The point at $n=k$ never plays a role to determine the mode. The joint modes
are at $0$ and the location of the left mountain peak, say $m_{\rm left}$. The
value of $m_{\rm left}$ depends on $\lambda$ but is always at least $k+2$. No
explicit formula is yet known for $m_{\rm left}$ when a double mode is
attained with the point at $n=0$, although the values can be tabulated for all
$k\in[15,41]$.

  4. 4.

When the value of $\lambda$ increases, the unique mode jumps from $0$ to
$m_{\rm left}$, i.e. an increase of more than one unit. This is the “first
double mode” (perhaps it is better to say the first mode jump) and is a
boundary between two domains in the parameter space.

  5. 5.

As the value of $\lambda$ increases further, the mode is determined by the
location of the left mountain peak. The value of the mode increases in unit
steps and there are double modes consisting of pairs of consecutive integers.
However, the mountain range develops a second (right) peak and its height
rises faster than that of the left peak and it catches up with the height of
the left peak.

  6. 6.

As the value of $\lambda$ increases further, a double mode is attained, where
the heights of the two mountain peaks are equal. This is displayed in the
fourth panel, where $\lambda\simeq 0.24159$. We denote the locations of the
two mountain peaks by $m_{\rm left}$ and $m_{\rm right}$, respectively. Very
significant: the value of $m_{\rm left}$ is larger than that in the previous
panel, i.e. the value of $m_{\rm left}$ increases with $\lambda$. Observe also
that there is a local minimum between the two mountain peaks, i.e. they are
separated by more than one unit and are not consecutive integers.

  7. 7.

When the value of $\lambda$ increases, the unique mode jumps from $m_{\rm
left}$ to $m_{\rm right}$, which is an increase of more than one unit. We
refer to this as the “second mode jump” (by which is meant an increase of the
mode by more than one unit) and is a boundary between two domains in the
parameter space.

  8. 8.

As the value of $\lambda$ increases further, the height of the right mountain
peak exceeds that of the left. The location of the right peak shifts
rightwards and its height increases. The mode is determined solely by the
location of the right mountain peak and increases in unit steps and takes all
integer values from its value at the second mode jump upwards.

  9. 9.

For a discrete (denumerably infinite) set of values of $\lambda$, the right
mountain peak has a flat top and the distribution is bimodal. The joint modes
consist of consecutive integers. This is displayed in the fifth panel, where
$\lambda\simeq 0.3039$. There is a double mode at $55$ and $56$.

###  3.5 Case $k\geq 42$

We select $k=50$ as a representative example. The scaled pmf of the Poisson
distribution of order $50$ is plotted in Fig. 6 for selected values of
$\lambda$, increasing from top to bottom.

  1. 1.

In the top panel, $\lambda=0.04$. This value is so small that the mountain
range region ($n>k$) is monotonically decreasing. The increasing sequence for
$n=1,\dots,k$ is visible, as is the local maximum at $n=k$, but it is not a
global maximum. The unique mode is zero.

  2. 2.

As the value of $\lambda$ increases, the mountain range develops a left peak
and its height rises and catches up with that of the point at $n=k$. This is
displayed in the second panel, where $\lambda\simeq 0.07822$. Note that both
heights are less than $1$, hence they are not modes.

  3. 3.

As the value of $\lambda$ increases further, the mountain range develops a
right peak and its height rises and catches up with that of the left peak.
This is displayed in the third panel, where $\lambda\simeq 0.098$. Note that
both heights are less than $1$, hence they are not modes.

  4. 4.

As the value of $\lambda$ increases further, the height of the right mountain
peak rises faster than that of the left peak and it equals $1$ before the left
peak or the point at $n=k$ does so. This is displayed in the fourth panel,
where $\lambda\simeq 0.10194$. Hence the first double mode is not at $0$ and
$k$. The point at $n=k$ and the left mountain peak never play a role to
determine the mode. The joint modes are at $0$ and the location of the right
mountain peak, say $m_{\rm right}$. The value of $m_{\rm right}$ depends on
$\lambda$ but is always at least $k+2$. No explicit formula is yet known for
$m_{\rm right}$ when a double mode is attained with the point at $n=0$,
although the values can be tabulated for all desired values $k\geq 42$.

  5. 5.

When the value of $\lambda$ increases, the unique mode jumps from $0$ to
$m_{\rm right}$, i.e. an increase of more than one unit. This is the “first
double mode” (perhaps it is better to say the first mode jump) and is a
boundary between two domains in the parameter space.

  6. 6.

As the value of $\lambda$ increases further, the mode is determined solely by
the location of the right mountain peak. The location of the right peak shifts
rightwards and its height increases. The mode increases in unit steps and
takes all integer values from its value at the first mode jump upwards. There
is only one mode jump, for $k\geq 42$.

  7. 7.

For a discrete (denumerably infinite) set of values of $\lambda$, the right
mountain peak has a flat top and the distribution is bimodal. The joint modes
consist of consecutive integers. This is displayed in the fifth panel, where
$\lambda\simeq 0.105$. There is a double mode at $116$ and $117$.

###  3.6 Parameter map

The shape of the pmf of the Poisson distribution of order $k$ can be
partitioned into three sections: (i) the single point at $n=0$, (ii) the
increasing sequence ($n\in[1,k]$), and (iii) the mountain range ($n>k$). It is
by no means obvious from the formal sum in eq. (1.1) that such a partition
exists. The increasing sequence is peculiar to the Poisson distribution of
order $k$ and does not exist for the standard Poisson distribution ($k=1$). It
consists of terms where the number of summands equals $n$, because
$n\in[1,k]$. The mountain range consists of terms where the number of summands
is strictly less than $n$, because $n>k$. The mountain range is similar in
character to the scaled pmf of the standard Poisson distribution. For
sufficiently small $\lambda>0$, it decreases monotonically, but for larger
values of $\lambda$ it exhibits a peak. It actually exhibits a maximum of two
peaks. The reason for this is not known: why not a single peak, and why not
more than two peaks?

The mode structure of the scaled pmf of the Poisson distribution of order $k$
is determined by four parameters: (i) the single point at $n=0$, whose height
is always $1$, (ii) the single point at $n=k$, whose height depends on
$\lambda$, (iii) the left mountain peak $m_{\rm left}$, and (iv) the right
mountain peak $m_{\rm right}$. For $k=2$ and $k=3$, there is only a single
peak in the region $n>k$. Both the location and the height of the mountain
peaks depend on $\lambda$. The evidence presented in this note suggests that
the Poisson distribution of order $k$ does not have three or more joint modes.
The evidence is admittedly numerical, and cannot claim to be exhaustive, hence
the existence of triple, etc. modes remains open.

This is the parameter map:

  1. 1.

For $k\in[2,3]$, the following happens as the value of $\lambda$ increases
continuously from $0$.

    1. (a)

For sufficiently small $\lambda>0$, the mode is zero.

    2. (b)

As $\lambda$ increases, the height of the point at $n=k$ reaches $1$ and there
is a double mode at $0$ and $k$.

    3. (c)

When this domain boundary is crossed, the mode jumps from $0$ to $k$. This is
the first mode jump.

    4. (d)

As $\lambda$ increases, the height of the single peak in the mountain range
catches up to the height of the point at $n=k$ and there is a double mode at
$k$ and $k+2$.

    5. (e)

When this domain boundary is crossed, the mode jumps from $k$ to $k+2$. This
is the second mode jump.

    6. (f)

As $\lambda$ increases further, the mode is determined by the single peak in
the mountain range. The location of the single peak shifts rightwards and its
height increases. The mode increases in unit steps and takes all integer
values $\geq k+2$. There is a denumerable infinity of double modes, consisting
of pairs of consecutive integers.

  2. 2.

For $k\in[4,14]$, the following happens as the value of $\lambda$ increases
continuously from $0$.

    1. (a)

For sufficiently small $\lambda>0$, the mode is zero.

    2. (b)

As $\lambda$ increases, the height of the point at $n=k$ reaches $1$ and there
is a double mode at $0$ and $k$.

    3. (c)

When this domain boundary is crossed, the mode jumps from $0$ to $k$. This is
the first mode jump.

    4. (d)

As $\lambda$ increases, the height of the left mountain peak catches up to the
height of the point at $n=k$ and there is a double mode at $k$ and $m_{\rm
left}$. The value of $m_{\rm left}$ is not less than $k+2$.

    5. (e)

When this domain boundary is crossed, the mode jumps from $k$ to $m_{\rm
left}$. This is the second mode jump.

    6. (f)

As $\lambda$ increases further, the mode is determined by the left mountain
peak, which shifts rightwards and its height increases. The mode increases in
unit steps and there are some double modes, consisting of pairs of consecutive
integers.

    7. (g)

As $\lambda$ increases further, the height of the right mountain peak catches
up to the height of the left mountain peak and there is a double mode at
$m_{\rm left}$ and $m_{\rm right}$. The values of $m_{\rm left}$ and $m_{\rm
right}$ are never consecutive integers.

    8. (h)

When this domain boundary is crossed, the mode jumps from $m_{\rm left}$ to
$m_{\rm right}$. This is the third mode jump.

    9. (i)

As $\lambda$ increases further, the mode is determined by the right mountain
peak, which shifts rightwards and its height increases. The mode increases in
unit steps and takes all values from the third mode jump upwards. There is a
denumerable infinity of double modes, consisting of pairs of consecutive
integers.

  3. 3.

For $k\in[15,41]$, the following happens as the value of $\lambda$ increases
continuously from $0$.

    1. (a)

For sufficiently small $\lambda>0$, the mode is zero.

    2. (b)

As $\lambda$ increases, the height of the left mountain peak catches up to the
height of the point at $n=k$, but their heights are less than $1$. The mode is
zero.

    3. (c)

As $\lambda$ increases further, the height of the left mountain peak reaches
$1$ and there is a double mode at $0$ and $m_{\rm left}$. The value of $m_{\rm
left}$ is not less than $k+2$. The point at $n=k$ plays no role to determine
the mode structure.

    4. (d)

When this domain boundary is crossed, the mode jumps from $0$ to $m_{\rm
left}$. This is the first mode jump.

    5. (e)

As $\lambda$ increases further, the mode is determined by the left mountain
peak, which shifts rightwards and its height increases. The mode increases in
unit steps and there are some double modes, consisting of pairs of consecutive
integers.

    6. (f)

As $\lambda$ increases further, the height of the right mountain peak catches
up to the height of the left mountain peak and there is a double mode $m_{\rm
left}$ and $m_{\rm right}$. The values of $m_{\rm left}$ and $m_{\rm right}$
are never consecutive integers.

    7. (g)

When this domain boundary is crossed, the mode jumps from $m_{\rm left}$ and
$m_{\rm right}$. This is the second mode jump.

    8. (h)

As $\lambda$ increases further, the mode is determined by the right mountain
peak, which shifts rightwards and its height increases. The mode increases in
unit steps and takes all values from the second mode jump upwards. There is a
denumerable infinity of double modes, consisting of pairs of consecutive
integers.

  4. 4.

For $k\geq 42$, the following happens as the value of $\lambda$ increases
continuously from $0$.

    1. (a)

For sufficiently small $\lambda>0$, the mode is zero.

    2. (b)

As $\lambda$ increases, the height of the left mountain peak catches up to the
height of the point at $n=k$, but their heights are less than $1$. The mode is
zero.

    3. (c)

As $\lambda$ increases further, the height of the right mountain peak catches
up to the height of the left mountain peak but their heights are less than
$1$. The mode is zero.

    4. (d)

As $\lambda$ increases further, the height of the right mountain peak reaches
$1$ and there is a double mode at $0$ and $m_{\rm right}$. The value of
$m_{\rm right}$ is not less than $k+2$. The point at $n=k$ and the left
mountain peak play no role to determine the mode structure.

    5. (e)

When this domain boundary is crossed, the mode jumps from $0$ to $m_{\rm
right}$. This is the first (and only) mode jump.

    6. (f)

As $\lambda$ increases further, the mode is determined by the right mountain
peak, which shifts rightwards and its height increases. The mode increases in
unit steps and takes all values from the first mode jump upwards. There is a
denumerable infinity of double modes, consisting of pairs of consecutive
integers.

###  3.7 Excluded values

We can now explain the “excluded values” tabulated in [2], i.e. integers which
cannot be modes of the Poisson distribution of order $k$. The values tabulated
in [2] are the integers which are skipped in the first, second or third mode
jumps, and we can now explain the cause of those mode jumps. A mode jump
greater than unity occurs only when the controlling parameter of the mode
changes, e.g. from $0$ to $k$ or from the left peak to the right peak, etc.
Only values $k\in[4,14]$ display three mode jumps, and for $k\geq 42$ there is
only one mode jump. But for all $k\geq 2$ there is always a mode jump. As
noted in [2], the integers $1$ and $k+1$ are never modes of the Poisson
distribution of order $k\geq 2$.

##  4 Median

It was proved in [1] that the median is zero if and only if $\lambda\leq(\ln
2)/k$. Numerical studies reported in [1] also yielded the following expression
for the median for $\lambda\geq 1$. If $n\in\mathbb{N}$ and $n\geq\kappa$, set
$\lambda=n/\kappa$ (so $\lambda\geq 1$). Then the median is given by (eq.
(3.1) in [1])

$\nu_{k}(n/\kappa)=n-\biggl{\lfloor}\frac{k+4}{8}\biggr{\rfloor}\,.$ (4.1)

Note that $\frac{\kappa}{k}\,\ln 2=\frac{1}{2}(k+1)\ln
2>\lfloor(k+4)/8\rfloor$ for all $k\geq 1$. Hence we conjecture the following
bounds for the value of the median in the intermediate zone
$\lambda\in(\frac{1}{k}\,\ln 2,\,1)$.

######  Conjecture 4.1.

For fixed $k\geq 2$ and $\lambda\in(\frac{1}{k}\,\ln 2,\,1)$, we claim the
lower bound for the median is

$\nu_{k}(\lambda)\geq\max\bigl{\\{}0,\lfloor\kappa\lambda\rfloor-{\textstyle\frac{1}{2}}(k+1)\ln
2\bigr{\\}}\,.$ (4.2)

For an upper bound, we know from the numerical studies reported in [1] it is
possible for the value of the median to exceed the mean, because
$\lfloor(k+4)/8\rfloor=0$ if $k<4$. We propose the following upper bound for
the median, where $c_{0}=1$ if $k\leq 3$ and $c_{0}=0$ if $k\geq 4$.

$\nu_{k}(\lambda)\leq\lfloor\kappa\lambda\rfloor+c_{0}\,.$ (4.3)

Monte Carlo scans using values $2\leq k\leq 4\times 10^{4}$ and
$\lambda\in(\frac{1}{k}\,\ln 2,\,1)$ found no violations of the bounds in eqs.
(4.2) or (4.3). Note the following:

  1. 1.

From eq. (4.1), it is tempting to conjecture that a sharper upper bound is
$\nu_{k}(\lambda)\leq\lfloor\kappa\lambda\rfloor-\lfloor(k+4)/8\rfloor$ or
possibly $\nu_{k}(\lambda)\leq
1+\lfloor\kappa\lambda\rfloor-\lfloor(k+4)/8\rfloor$, but they are both false
for $\lambda<1$. An example is $k=514$ and $\lambda\simeq 0.0031619$, whence
$\kappa\lambda\simeq 418.4998$ and
$\lfloor\kappa\lambda\rfloor-\lfloor(k+4)/8\rfloor=418-64=354$, but the median
value is $367$.

  2. 2.

A graph of the median is plotted in Fig. 7, for $k=10$ and
$0<\kappa\lambda\leq 40$, whence $\kappa=55$ and $0<\lambda\leq
40/55=0.7272\dots$ The value of the median is plotted as the solid line. The
lower bound from eq. (4.2) is plotted as the dotted line. The upper bound from
eq. (4.3) is plotted as the dashed line. Observe that the upper bound equals
the median at several places, indicating that it is a sharp upper bound.
Observe also that the lower bound is essentially a straight line, and does not
capture the curvature as the value of the median approaches zero for small
values of $\kappa\lambda$. This indicates that the expression for the lower
bound can be improved.

  3. 3.

The graph of the median in Fig. 7 looks hyperbolic, with a straight line
asymptote for large values of $\kappa\lambda$ (see eq. (4.1) for $\lambda\geq
1$). This suggests a better parameterization for the value of the median is
possible in the interval $\lambda\in(\frac{1}{k}\ln 2,\,1)$. The matter is
left for future work.

Given the above, we propose the following exact results and bounds for the
median, for all $k\geq 2$ and $\lambda>0$. In all cases, we fix $k\geq 2$ and
the median is denoted by $\nu_{k}(\lambda)$.

  1. 1.

For $\lambda\in(0,\frac{1}{k}\,\ln 2]$ the median is zero:
$\nu_{k}(\lambda)=0$. This was proved in [1].

  2. 2.

For $\lambda\in(\frac{1}{k}\,\ln 2,\,1)$, the median is bounded via eqs. (4.2)
and (4.3), where $c_{0}=1$ if $k\leq 3$ and $c_{0}=0$ if $k\geq 4$.

$\max\bigl{\\{}0,\lfloor\kappa\lambda\rfloor-{\textstyle\frac{1}{2}}(k+1)\ln
2\bigr{\\}}\leq\nu_{k}(\lambda)\leq\lfloor\kappa\lambda\rfloor+c_{0}\,.$ (4.4)

  3. 3.

For $\lambda\geq 1$, the median is given by eqs. (3.1), (3.3) and (3.4) in [1]
as follows. Let $n\in\mathbb{N}$ and $n\geq\kappa$, then for
$\kappa\lambda\in(\alpha_{k,n-1},\alpha_{k,n}]$ the median equals
$\nu_{k}(n/\kappa)$ as follows.

$\displaystyle\nu_{k}(n/\kappa)$
$\displaystyle=n-\biggl{\lfloor}\frac{k+4}{8}\biggr{\rfloor}\,,$ (4.5a)
$\displaystyle\alpha_{k,n}$
$\displaystyle=n+\textrm{frac}\Bigl{(}\frac{k+4}{8}\Bigr{)}+\frac{k}{8(2k+1)}+A_{k,n}\,,$
(4.5b) $\displaystyle A_{k,n}$
$\displaystyle=\biggl{(}\frac{3\kappa}{349}+\frac{13}{1000}\biggr{)}\frac{1}{n}+\frac{13}{1500}\biggl{(}\biggl{\lfloor}\frac{k+4}{8}\biggr{\rfloor}-3\biggr{)}\,\frac{\kappa}{n^{2}}+\cdots$
(4.5c)

The expression for $A_{k,n}$ is approximate, but the expressions for
$\nu_{k}(n/\kappa)$ and $\alpha_{k,n}$ are otherwise conjectured to be exact.

##  5 Mode

Recall $\hat{m}_{k}$ is the value of the first double mode of the Poisson
distribution of order $k$, i.e. the distribution is bimodal, with modes at $0$
and $\hat{m}_{k}$ and $\hat{\lambda}_{k}$ is the corresponding value of
$\lambda$. It was shown in [2] that $\hat{\lambda}_{k}$ is a strictly
decreasing function of $k$ and also that $\hat{m}_{k}\geq k$.

The following asymptotic expression for the mode was conjectured in [1]. Let
$n\in\mathbb{N}$ and $n\geq 2\kappa$ and set $\lambda=n/\kappa$ (so
$\lambda\geq 2$). Then the mode is given by (eq. (4.1) in [1])

$m_{k}(n/\kappa)=n-\biggl{\lfloor}\frac{3k+5}{8}\biggr{\rfloor}\,.$ (5.1)

Theorem 2.1 in [8] states the following upper and lower bounds for the mode.

$\bigl{\lfloor}\kappa\lambda\bigr{\rfloor}-\kappa+1-\delta_{k,1}\leq
m_{k}(\lambda)\leq\bigl{\lfloor}\kappa\lambda\bigr{\rfloor}\,.$ (5.2)

We conjecture an improved lower bound for the mode.

######  Conjecture 5.1.

For fixed $k\geq 2$ and $\lambda\in(\hat{\lambda}_{k},\,2)$ (so the value of
the mode is nonzero), we propose the following as an improved lower bound for
the mode

$m_{k}(\lambda)\geq\lfloor\kappa\lambda\rfloor-k\,.$ (5.3)

The right-hand side is nonnegative and evidence will be presented below that
it is a sharp lower bound.

Monte Carlo scans using values $2\leq k\leq 4\times 10^{4}$ and $0<\lambda<2$
found no violations of eq. (5.3). (Only cases where the value of the mode was
positive were included in the scan, to satisfy the requirements of Conjecture
5.1.) Note the following:

  1. 1.

From eq. (5.1), it is tempting to conjecture that a sharper upper bound is
$m_{k}(\lambda)\leq\lfloor\kappa\lambda\rfloor-\lfloor(3k+5)/8\rfloor$ or
possibly $m_{k}(\lambda)\leq
1+\lfloor\kappa\lambda\rfloor-\lfloor(3k+5)/8\rfloor$, but they are both false
for $\lambda<2$. An example is $k=44$ and $\lambda\simeq 0.114198$, whence
$\kappa\lambda\simeq 113.056$ and
$\lfloor\kappa\lambda\rfloor-\lfloor(3k+5)/8\rfloor=113-17=96$, but the mode
value is $98$.

  2. 2.

Note however that Monte Carlo scans have thus far failed to find an example
where $m_{k}(\lambda)>0$ and $m_{k}(\lambda)\geq
3+\lfloor\kappa\lambda\rfloor-\lfloor(3k+5)/8\rfloor$.

  3. 3.

A graph of the mode is plotted in Fig. 8, for $k=10$ and $0<\kappa\lambda\leq
40$, whence $\kappa=55$ and $0<\lambda\leq 40/55=0.7272\dots$ The value of the
mode is plotted as the solid line. The lower bound from eq. (5.3) is plotted
as the dotted line. The upper bound is $\lfloor\kappa\lambda\rfloor$ (from
[8]) and is plotted as the dashed line.

  4. 4.

Unlike the median, the value of the mode does not always increase in unit
steps. Observe in Fig. 7 that the value of the mode jumps from $0$ to $10$ and
then from $10$ to $17$ and then from $20$ to $23$. These mode jumps were
explained in Sec. 3.

  5. 5.

In Fig. 7, it looks as if the lower bound equals the mode at
$\kappa\lambda=20$, but it does not. The value of the mode jumps from $10$ to
$17$ at $\kappa\lambda\simeq 19.91$, while the lower bound increases from $9$
to $10$ at $\kappa\lambda=20$.

  6. 6.

However, the case $k=2$ demonstrates that the expression for the lower bound
in eq. (5.3) is a sharp lower bound. A graph of the mode is plotted in Fig. 9,
for $k=2$ and $0<\kappa\lambda\leq 8$, whence $\kappa=3$ and $0<\lambda\leq
8/3=2.666\dots$ The value of the mode is plotted as the solid line. The upper
and lower bounds are plotted as the dashed and dotted lines, respectively.
Both the mode and lower bound equal $2$ at $\kappa\lambda=4$. For
$\kappa\lambda=4.0238$ the value of the mode jumps to $4$ but the lower bound
remains at $2$. Hence the mode and lower bound both equal $2$ in a nonempty
subset of the interval $\kappa\lambda\in[4,4.0238)$.

Given the above, we propose the following exact results and bounds for the
mode, for all $k\geq 2$ and $\lambda>0$. In all cases, we fix $k\geq 2$ and
the mode is denoted by $m_{k}(\lambda)$.

  1. 1.

For $\lambda\in(0,\hat{\lambda}_{k})$ the mode is zero: $m_{k}(\lambda)=0$, by
the definition of $\hat{\lambda}_{k}$.

  2. 2.

For $\lambda=\hat{\lambda}_{k}$, there is a double mode at $0$ and
$\hat{m}_{k}$ (this is the definition of $\hat{m}_{k}$).

  3. 3.

For $\lambda\in(\hat{\lambda}_{k},\,2)$, the mode is bounded as follows

$\lfloor\kappa\lambda\rfloor-k\leq
m_{k}(\lambda)\leq\lfloor\kappa\lambda\rfloor\,.$ (5.4)

The lower bound is based on numerical studies in this note. The upper bound is
from [8].

  4. 4.

For $\lambda\geq 2$, the mode is given by eqs. (4.1), (4.5) and (4.6) in [1]
as follows. Let $n\in\mathbb{N}$ and $n\geq 2\kappa$, then for
$\kappa\lambda\in(\beta_{k,n-1},\beta_{k,n}]$ the mode equals
$m_{k}(n/\kappa)$ as follows.

$\displaystyle m_{k}(n/\kappa)$
$\displaystyle=n-\biggl{\lfloor}\frac{3k+5}{8}\biggr{\rfloor}\,,$ (5.5a)
$\displaystyle\beta_{k,n}$
$\displaystyle=n+\textrm{frac}\Bigl{(}\frac{3k+5}{8}\Bigr{)}+\frac{k-1}{8(2k+1)}+B_{k,n}\,,$
(5.5b) $\displaystyle B_{k,n}$
$\displaystyle=\biggl{(}\frac{\kappa}{16+\frac{8}{9}}-\frac{1}{13+\frac{2}{3}}\biggr{)}\frac{1}{n}+\biggl{\lfloor}\frac{3k+5}{8}\biggr{\rfloor}\,\frac{3\kappa}{50n^{2}}+\cdots$
(5.5c)

The expression for $B_{k,n}$ is approximate, but the expressions for
$m_{k}(n/\kappa)$ and $\beta_{k,n}$ are otherwise conjectured to be exact.

  5. 5.

Note that no general formula is yet known for $\hat{m}_{k}$ or
$\hat{\lambda}_{k}$, although specific cases have been solved.

It was stated in [1] that the asymptotic formula for the mode is applicable
for values $\lambda\geq 2$ (see eq. (5.1)). Just out of curiosity, let us plot
the pmf for $\lambda=2$ for the values $k=3,10,20,50$ employed in Sec. 3. To
plot them on the same scale, we scale all the curves to a peak height of $1$
and on the horizontal axis we plot the value of $n/\kappa\in[0,6]$. The curves
are displayed in Fig. 10. The mean is $\kappa\lambda=2\kappa$, hence the mean
is at $n/\kappa=2$ for all the curves. Observe that all the curves are smooth.
Issues of a local maximum at $n=k$, or a left peak, have faded out of the
picture. The smoothness improves as $k$ increases, which is expected. The
peaks get narrower as $k$ increases. We can explain this as follows. Recall
the variance is $\sigma_{k}^{2}(\lambda)=\frac{1}{6}k(k+1)(2k+1)\lambda$,
hence in Fig. 10 the scaled peak width is (scaled standard deviation)

$\begin{split}\frac{\sigma_{k}(\lambda)}{\kappa}&=\sqrt{\frac{\frac{1}{6}k(k+1)(2k+1)\lambda}{\frac{1}{4}k^{2}(k+1)^{2}}}\\\
&=\sqrt{\frac{2\lambda}{3}}\,\sqrt{\frac{2k+1}{k(k+1)}}\\\
&=\sqrt{\frac{2\lambda}{3}}\,\sqrt{\frac{1}{k}+\frac{1}{k+1}}\;.\end{split}$
(5.6)

Hence the scaled peak width decreases as the value of $k$ increases. We fixed
$\lambda=2$ in Fig. 10. Observe also that the peaks shift to the right as $k$
increases. This is consistent with eq. (5.1). Set
$n=\kappa\lambda=\mu_{k}(\lambda)$ and divide by $\kappa$ to obtain

$\begin{split}\frac{\mu_{k}(\lambda)-m_{k}(\lambda)}{\kappa}&=\frac{\lfloor(3k+5)/8\rfloor}{\frac{1}{2}k(k+1)}\\\
&\leq\frac{(3k+5)/8}{\frac{1}{2}k(k+1)}\\\
&=\frac{1}{4}\,\biggl{(}\frac{5}{k}-\frac{2}{k+1}\biggr{)}\,.\end{split}$
(5.7)

The scaled difference between the mean and the mode decreases as the value of
$k$ increases.

##  6 Conclusion

The major goal of this note was to characterize the structure of the
probability mass function (pmf) of the Poisson distribution of order $k$. The
pmf can be partitioned into a single point at $n=0$, an increasing sequence
for $n\in[1,k]$ and a mountain range for $n>k$. That structure is by no means
obvious from the formal definition of the pmf as a sum over terms with
factorial denominators. The mode structure of the pmf was quantified. The
locations of and the reasons for the mode jumps (where the mode increases by
more than one unit) were established. The “excluded values” tabulated in [2],
i.e. integers which cannot be modes of the Poisson distribution of order $k$,
were explained as those integers which are skipped in the mode jumps. It was
demonstrated that for all $k\geq 2$, the Poisson distribution of order $k$ has
a denumerably infinite set of double modes, consisting of pairs of consecutive
integers. Numerical evidence was presented that the Poisson distribution of
order $k$ does not have three or more joint modes. The opportunity was also
taken to publish improvements to various inequalities (sharper bounds, etc.)
and also to present new conjectured upper and lower bounds for the median and
the mode of the Poisson distribution of order $k$.

## References

  * [1] S.R. Mane, “Asymptotic results for the Poisson distribution of order $k$” arXiv:2309.05190 [math.PR] (2023).
  * [2] S.R. Mane, “First double mode of the Poisson distribution of order $k$” arXiv:2309.09278 [math.PR] (2023).
  * [3] A.N. Philippou, C. Georghiou and C. Philippou, “A generalized geometric distribution and some of its properties” Stat. Prob. Lett. 1, 171–175 (1983).
  * [4] R. M. Adelson, “Compound Poisson Distributions” Operational Research Quarterly 17, 73–75 (1966).
  * [5] A.N. Philippou, “Poisson and compound Poisson distributions of order $k$ and some of their properties” Journal of Soviet Mathematics 27, 3294–3297 (1984).
  * [6] J.A. Adell and P. Jodrá, “The median of the Poisson distribution” Metrika 61, 337–346 (2005).
  * [7] Y. Kwon and A.N. Philippou, “The Modes of the Poisson Distribution of Order 3 and 4” Entropy 25, 699 (2023).
  * [8] C. Georghiou, A.N. Philippou and A. Saghafi, “On the Modes of the Poisson Distribution of Order $k$” Fibonacci Quarterly 51, 44–48 (2013).
  * [9] A.N. Philippou, “a note on the modes of the poisson distribution of order $k$” Fibonacci Quarterly 52, 203–205 (2014).

Figure 1:  Plot of the scaled pmf $h_{k}(n;\lambda)$ of the standard Poisson
distribution ($k=1$) for various values of the rate parameter $\lambda$.

Figure 2:  Plot of the scaled pmf $h_{k}(n;\lambda)$ for the Poisson
distribution of order $50$ with rate parameter $\lambda=0.10194$.

Figure 3:  Plot of the scaled pmf $h_{k}(n;\lambda)$ for the Poisson
distribution of order $3$ for selected values of the rate parameter $\lambda$.

Figure 4:  Plot of the scaled pmf $h_{k}(n;\lambda)$ for the Poisson
distribution of order $10$ for selected values of the rate parameter
$\lambda$.

Figure 5:  Plot of the scaled pmf $h_{k}(n;\lambda)$ for the Poisson
distribution of order $20$ for selected values of the rate parameter
$\lambda$.

Figure 6:  Plot of the scaled pmf $h_{k}(n;\lambda)$ for the Poisson
distribution of order $50$ for selected values of the rate parameter
$\lambda$.

Figure 7:  Plot of the median (solid line) for the Poisson distribution of
order $10$, with the upper bound (dashed) and lower bound (dotted).

Figure 8:  Plot of the mode (solid line) for the Poisson distribution of order
$10$, with the upper bound (dashed) and lower bound (dotted).

Figure 9:  Plot of the mode (solid line) for the Poisson distribution of order
$2$, with the upper bound (dashed) and lower bound (dotted).

Figure 10:  Plot of the pmf for the Poisson distribution of order $k$ for
$\lambda=2$ and $k=3,10,20,50$. For each value of $k$, the pmf was scaled to a
peak height of $1$.


# Zinc dialkyldithiophosphates adsorption and dissociation on ferrous
substrates:
an ab initio study

Francesca Benini  Paolo Restuccia  Maria Clelia Righi<EMAIL_ADDRESS>Physics and Astronomy Department, University of Bologna,
Viale Berti Pichat 6/2, 40137, Bologna (Italy)

###### Abstract

Zinc dialkyldithiophosphates (ZDDPs) have been commonly used as anti-wear
additives in the automotive industry for the past 80 years. The morphology,
composition and structure of the ZDDPs phosphate-based tribofilm, which is
essential for its lubricant functioning, have been widely studied
experimentally. However, despite their widespread use, a general agreement on
their primary functioning mechanism is still lacking. The morphology and
composition of the ZDDPs phosphate-based tribofilm have been widely studied
experimentally, but the formation process and the relevant driving forces are
still largely debated. In particular, it is unclear whether the stress-induced
molecular dissociation occurs in the bulk oil or on the substrate. In this
work, we employ ab initio density-functional theory simulations to compare
ZDDP fragmentation in vacuum and over a reactive substrate, considering the
effects of surface oxidation on the dissociation path. Our results show that
the molecular dissociation is highly endothermic in the absence of a
supporting substrate, while in the presence of an iron substrate it becomes
highly energetically favoured. Moreover, the presence of the substrate changes
the reaction path. At the same time, surface oxidation reduces the molecule-
substrate interaction. These findings provide valuable insights into the early
stages of the formation of phosphate-based tribofilms.

##  I Introduction

After more than 80 years since their discovery, zinc dialkyldithiophosphates
(ZDDPs) are still the most used anti-wear additives in the automotive industry
[1, 2, 3]. Indeed, this class of materials is known for their capability to
prevent wear, especially on steel [4, 5, 6], thanks to the formation of a
glassy phosphate-based tribofilm [7, 1].

Many details on the ZDDP tribofilm chemical composition and morphology are
known thanks to the great effort that has been put experimentally into its
characterization[8, 9, 7, 10, 11, 12, 13, 14]: the film, typically around
80-150 nm thick and composed of flat-topped pads [10], consists of a
superficial layer of zinc and iron polyphosphate and/or phosphate; closer to
the steel surface it presents a sulfur-rich layer [1, 8, 13]; finally,
although believed amorphous for years, it was recently found that ZDDP
tribofilms hide a nanocrystalline structure, promoted by prolonged rubbing,
which improves their durability. [11]

Concerning the tribofilm formation process, there is a common agreement
regarding its main steps [11, 14]: first, surface adsorption and the initial
formation of iron disulphide take place; secondly, the alkyl group transfers
from oxygen to sulfur, converting dialkylthiophosphate to dithionyl phosphate;
thirdly, intermolecular reactions enable phosphates to polymerize to build
polyphosphate chains, the building blocks for the tribofilm; and finally, the
amorphous film forms, which then becomes crystalline within a depolymerization
process enhanced by prolonged rubbing at elevated temperatures [11].

However, despite such general consensus, there are still matters of debate,
for instance concerning what specifically drives ZDDP dissociation and the
subsequent film formation, as well as the exact fragmentation path followed by
the molecule and the role played by the substrate in such process [15, 12,
13]. Concerning the latter, some authors proposed that the substrate hardness
might affect tribofilm quality [4, 16, 17]. Indeed, the appearance of a
protective film in tribological conditions is a common phenomenon, and many
factors might affect its formation. This is why a deep understanding of the
different roles played by various physical quantities, such as compressive and
shear stress, temperature, and the interaction with the substrate, is crucial
to design new lubricant additives [18] and unveil the functionality of
commercial ones [19]. These contributions are difficult to isolate
experimentally in actual tribological conditions since they are all at play
simultaneously at the nano-asperity contact.

Among all of these quantities, it has been suggested that shear stress within
the base oil might be the driving force for ZDDP tribofilm formation. This
hypothesis, described through the Stress-Promoted Thermal Activation model [7,
20, 8, 12], has been strengthened by experimental evidence of ZDDP tribofilm
formation on a variety of substrates, ranging from metals [21, 22] to ceramics
[23, 24, 25], DLC coatings [26] and silicon [4], suggesting that the molecule-
substrate chemical interaction plays a minor role in the phenomenon. However,
recent studies emphasized the role of reactive substrates, like the ones
exposed during tribological experiments, as the primary driving force in the
dissociation and tribofilm formation of commercial lubricant additives [27,
28, 18].

Finally, despite the widespread use of ZDDP additives, there is a scarcity of
theoretical atomistic simulations that can provide insights into the ZDDP
tribofilm formation process. In particular, the fundamental steps of the
chemical reaction paths leading to the ZDDP dissociation, inaccessible from
experiments, and the possible role of the molecule-substrate interaction are
still missing. Indeed, none of the limited numbers of computational works
concerning ZDDP takes into account the presence of a substrate [29, 30, 15,
31, 32]. In this work, we employ density-functional theory (DFT) to
investigate ZDDP fragmentation on iron and oxygen-passivated iron substrates
relevant to many technological applications. We study the dissociation path
followed by the ZDDP molecule when bonds are stretched in vacuum to model the
situation in which stresses arise in chemically non-reactive base oils and
compare it with molecule fragmentation on ferrous substrates. In this way, we
identified the most favorable dissociation path and analysed the effect of
substrate passivation on ZDDP adsorption.

##  II Methods

All calculations were performed using spin-polarized Density Functional Theory
(DFT) as implemented in version 6.8 of the Quantum Espresso suite [33, 34,
35]. The exchange-correlation functional was described using the generalized
gradient approximation (GGA) within the Perdew-Burke-Ernzerhof (PBE)
parametrization [36]. The kinetic energy cutoff for the wave functions was set
to 40 Ry, with a cutoff for charge densities of 320 Ry. A gaussian smearing
with a 0.001 Ry width was included to better describe occupations around the
Fermi level. The electronic configuration of atoms was described using
ultrasoft pseudopotentials in the RRKJ parametrization [37]. The default
energy and forces convergence threshold criteria (i.e., 10-4 Ry and 10-3
Ry/Bohr, respectively) were used to optimize the structural geometries. The
geometrical configuration and computational setup we used to describe the ZDDP
molecule, was already tested by our group in a previous work [38]. In
particular, the two lateral alkyl chains have been reduced to methyl groups to
minimize the computational cost. Here we study the interaction of the ZDDP
molecule with four different substrates: a clean Fe(110) four-layers thick
slab, and, by passivating it with oxygen atoms on both sides, oxidised Fe
slabs with 0.25 ML (Figure 1a), 0.5 ML (Figure 1b) and, 1 ML (Figure 1c)
O-coverages , respectively.

Figure 1: The oxidised substrates modelled in this study: 0.25ML (panel a),
0.5ML (panel b) and 1ML O-passivated Fe (panel c). From now on, Fe will be
represented in blue and O in red.

To obtain an accurate description of the most effective dissociative path and
of the role played by the substrate, the fragmentation of the molecule was
studied by bond stretching [39] comparing the results obtained in vacuum and
in the presence of the clean Fe substrate. To carry out this task in vacuum,
we placed the molecule vertically in a large supercell, and performed a
sequence of relax calculations keeping the central Zn atom fixed while
displacing firstly a P, then O, and finally a C atom of 0.2 Å respectively
along the $z$ direction at each step, reaching a total displacement of 5 Å.
The $z$ coordinates of the displaced atoms were kept fixed, and the optimized
coordinates of each step were used as the initial configuration for the next
step. An analogous procedure was performed for the case of the molecule on the
iron slab. We considered the optimized molecule adsorbed over the substrate as
a starting configuration. Then, we stretched the molecule by displacing first
P, then O and finally a C atom along the direction parallel to the surface,
namely the $x$ direction, keeping their $x$ coordinate fixed together with the
in-plane coordinates of the Zn atom. In this case, we reached a total
displacement of 3 Å, which was enough to see the complete fragmentation of the
molecule. We could then compute the reaction energies and force components
acting on the fixed atoms at each step along the $z$ ($x$) direction for the
stretching in the vacuum (Fe substrate) case. We expressed the obtained forces
as pressures dividing them by the in-plane area of the cell employed for the
molecule in vacuum (2.76 nm2) to allow for a rough estimation of the shear
stress that should be applied to obtain the corresponding molecular
deformation.

To study the adsorption of the ZDDP molecule on surfaces, we employed large 6
$\times$ 4 orthorhombic supercells containing 48 atoms per layer to avoid
undesired lateral interaction among the molecule replicas. Moreover, a
reasonable amount of void (at least 10 Å) was added to each cell to guarantee
a vertical distance among replicas. The resulting supercell dimensions were
17.10 Å $\times$ 16.12 Å $\times$ 30.78 Å. The latter was obtained after
optimizing a smaller 2 $\times$ 2 cell. Considering the large dimensions of
the supercell, we used a $\Gamma$ point sampling of the Brillouin zone. We
computed adsorption energies as follows:

$E_{ads}=E_{tot}-(E_{surf}+E_{adsorbate})$ (1)

where $E_{tot}$ is the energy of the optimized system composed of the
substrate and the adsorbed molecule/atom, while $E_{surf}$ ($E_{adsorbate}$)
is the energy of the isolated surface (adsorbate).

To study the ZDDP dissociation over iron substrates, we used the molecular
fragments identified by the bond stretching study. For each of them, we
identified the most favourable binding site and molecular orientation
employing Xsorb [40] and computing the dissociation energy $E_{frag}$ as
follows:

$E_{frag}=\sum_{n=1}^{N}E_{F_{n}}-(n\times E_{surf}+E_{mol})$ (2)

where $E_{F_{n}}$ represents the total energy of the most stable systems
composed of the surface and each fragment alone at a time, and $E_{surf}$
($E_{mol}$) represents the energy of the isolated surface (molecule). The
value of $n$ in the formula depends on how many fragments the molecule is cut
into. In our case, $n$ = 2, 3 and 4, respectively.

We further studied the ZDDP dissociative path looking at what happens at the
interface of two Fe (or 0.25 ML and 1 ML oxidised Fe) surfaces in the presence
of 1 GPa external load. To achieve this task, we realized uniaxial
compressions in which we reduced the vertical size of the cell and allowed the
ZDDP molecule to interact with one side of the slab and its own replica to
simulate the interface. This approach induced molecular dissociation and gave
us insights into the bond breaking that takes place at the interface under
load.

##  III Results and discussion

Here we report the results of our simulations: first, we analyse the effect of
iron oxidation on the adsorption of molecular ZDDP. Secondly, we compare ZDDP
dissociation in vacuum and in the presence of the clean substrate. We consider
that the dissociation in vacuum is comparable to the one within the base oil,
due to the shallow interaction between the lubricant additive and the
hydrocarbon molecule of the oil. We then look at the effect of substrate
oxidation on the dissociation energies, and finally consider how the
dissociation proceeds when the molecule is set at the interface between two
facing substrates.

Modeling bulk Fe oxides, like hematite, is extremely computationally
demanding, and it requires an exquisite fine-tuning of the computational
parameters, thus we modelled the effect of iron oxidation using iron slabs
passivated by oxygen atoms, as done in previous studies [19, 27, 38].

###  III.1 ZDDP adsorption

The first important step in the formation of ZDDP tribofilms is the adsorption
of the molecules on the substrate [1, 8, 14]. In addition to the clean iron
substrate we simulated three different degrees of surface O-passivation,
namely 0.25 ML, 0.5 ML and 1 ML, shown in Figure 1.

After fully optimizing the geometry of the adsorbed molecule/substrate system,
$E_{ads}$ was computed following Eq. 1. The optimized configurations,
alongside the corresponding adsorption energies, are reported in Fig. 2. For
all the considered systems, the adsorption energies are negative, indicating
that the molecular adsorption is energetically favourable. This is in
agreement with the experimental evidence that ZDDP undergoes chemisorption and
thermal decomposition [1, 41]. However, oxygen passivation significantly
reduces the interaction of the ZDDP molecule with Fe. In particular, $E_{ads}$
decreases by about one order of magnitude, from $-$1.03 eV for the bare iron
substrate to $-$0.12 and $-$0.11 eV for the half and fully passivated
substrates. The 0.25 ML passivated substrate maintains a mild reactivity, with
$E_{ads}$= $-$0.67 eV, as the ZDDP undergoes structural deformations leading
to bonds breaking and S interaction with iron atoms (Fig. 2b).

Figure 2: ZDDP adsorption on a) clean Fe and O passivated Fe with three
different O coverages: b) 0.25ML, c) 0.5ML and d) 1ML respectively. Adsorption
energies are also reported. From now on, P, Zn, C, H and S atoms are
represented in orange, grey, black, white and yellow, respectively.

###  III.2 Shear-induced dissociation of ZDDP in the absence and in the
presence of a substrate

Figure 3: Energy variation as a function of atom displacement along the $z$
direction (panel a) alongside force per unit area on P, C and O for the ZDDP
molecule stretching in vacuum (panel b). $F$ is the value of the force acting
on P, C and O whereas $A$ is the in-plane area of the cell. A positive
(negative) value of $F$ indicates a force with the same (opposite) direction
of the $z$ axis in the reference system.

The results of the bond-stretching procedure applied in vacuum are shown in
Fig. 3. In panel a) the energy variation as a function of the atom
displacement along the $z$ direction are displayed for all three pulling
atoms, whereas in panel b) the forces per unit area acting on P, C and O are
shown. In all the three considered cases, the energy increases during the atom
pulling until a plateau is reached at a fragments distance where the
interaction is negligible. Independently from the pulling atom, the detachment
of the S2P(OCH3)2 fragment has an associated energy cost of about 3.5 eV.
Alongside the energy, the force divided by the in-plane cell area acting on
the pulled atoms, provides a consistent picture with a maximum increase of the
resistant forces per unit area when bond breaking occurs, of about 1.2 GPa.

The geometry configurations along the stretching bond procedure are shown in
Fig. 4. Confirming what previously obtained in Ref. 38, we find that, whatever
atom is pulled, the most favourable fragmentation pattern involves breaking
the Zn-S bonds. Therefore, we found that in vacuum or inside the lubricating
oil, when shear stress is applied to the molecule, the Zn-S bonds are the less
energetically demanding to break, as also reported in the literature. [42, 7]

Figure 4: Configurations resulting from pulling, respectively, a P atom (panel
a), an O atom (panel b) and a C atom (panel c). The atom pulled is highlighted
with a dotted red circle and the displacement along $z$ is reported.

The results of the bond-stretching procedure on the iron substrate are
summarized in Fig. 5. Interestingly, the energy and the resistant stress in
the presence of a substrate are much lower than those in the absence,
indicating that the presence of the substrate makes molecular dissociation
much easier than in the oil. Indeed, the dissociated fragments are stabilized
by the substrate, as visible from the reaction energies that diminished
significantly to around $-6$ eV for both C and O pulling (Fig. 5a). In the
case of P, the reaction energy is significantly larger due to the shallower
interaction of the Zn-S central unit compared to the other two cases (as shown
in Fig. 6a).

Figure 5: Energy variation as a function of atom displacement along the $x$
direction (panel a) alongside force per unit area on P, C and O for the ZDDP
molecule stretching on an iron substrate (panel b). $F$ is the absolute value
of the force acting on P, C and O whereas $A$ is the in-plane area of the cell
employed for the molecule in vacuuum, so that values are comparable. A
positive (negative) value of $F$ indicates a force with the same (opposite)
direction of the $x$ axis in the reference system.

Not only the energies and forces are dramatically changed by the presence of
the substrate, but also the dissociation path is different. As shown in Fig.
6, the most favourable dissociation path involves the breaking of the P-S
bonds. Therefore, the P(OCH3)2 groups are released and chemisorbed to the iron
substrate, leaving the central Zn-S units and isolated S atoms adsorbed on
iron. These results are in agreement with the experimental observation of an
S-rich iron layer and phosphate layers within the tribofilm. A change in the
dissociative path induced by the presence of a reactive substrate was found
for MoDTC [19], another largely employed commercial lubricant additive. Our
results suggest that the substrate-mediated dissociation is much more probable
than the suggested shear-induced fragmentation within the oil [7] because the
energies and forces involved are significantly lower. Moreover, the reactive
substrate allowed to isolate the organophosphorus units from the central Zn-S
ones, thus explaining the evidence of separate phosphate and zinc sulfide
areas observed experimentally [1, 13].

Figure 6: ZDDP pulling on clean Fe. The atoms pulled, respectively P (panel
a), C (panel b) and O (panel c) are highlighted with a dotted red circle.
Atoms displacement is reported.

Starting from the hypothesis that the molecule may dissociate in bulk oil, we
first simulated the adsorption of the two fragments resulting from bond
stretching in vacuum over a substrate (Fig. 7). We performed a systematic
study of their adsorption by employing Xsorb [40], a computational tool
developed in our group, to perform a configurational study testing different
fragment orientations for each system. For the substrates, we take into
account the two opposite conditions of O-passivation, namely a bare Fe(110)
surface and a complete oxygen-passivated (i.e., 1 ML) Fe(110) substrate. We
found the optimal adsorption geometries reported in Fig. 7, together with the
corresponding $E_{frag}$.

Figure 7: ZDDP fragments adsorption following a two- (panel a), three- (panel
b) and four-fragments (panel c) dissociation path. In panel a) both adsorption
on 1 ML O-passivated Fe and clean Fe are reported, while panel b and c involve
bare iron only.

From this analysis shown in Fig. 7a, we found that the adsorption of the Zn-
terminated (F1) and the organophosphorus fragments (F2) over the O-passivated
substrate is endothermic, with an associated energy cost of $+$0.59 eV.
However, the same dissociation path on the bare iron substrate is
energetically favorable, with $E_{frag}$ equal to $-$3.86 eV. These results
confirm the effectiveness of O-passivation in hindering the ZDDP fragmentation
completely. We highlight that the F2 adsorbed over the clean surface caused
the S atoms (F2.a fragment) to further dissociate from the organophosphorus
group (F2.b fragment). We tested this hypothesis by considering a three-
fragments dissociation path (Fig. 7b), which resulted in a further $E_{frag}$
decrease to $-$5.37 eV. This result suggests that organophosphorus groups
prefer to stay isolated from the central units when adsorbed on iron, which is
consistent with previous results of similar compounds on iron 43. Complete
fragmentation of the organophosphorus groups from the ZDDP (Fig. 7c) resulted
in the lowest $E_{frag}$ equal to $-$8.61 eV. We gained useful knowledge by
studying the energeticis involved in dissociating ZDDP over a reactive
substrate. Specifically, we found that it is energetically favorable to
separate the P-terminated and Zn-S groups. This dissociation is the first step
in creating the zinc sulfide and phosphate layers seen in previous experiments
[1, 13].

###  III.3 Load-induced ZDDP dissociation

To verify the effects of load on molecular dissociation, we performed uniaxial
compressions of the molecule confined at an iron interface in the presence of
1 GPa load, as explained in the Methods section.

In Fig. 8, snapshots of the relaxation runs while reaching the target pressure
are shown for the different oxidation levels. In Fig. 8a and b, the snapshots
relative to the clean and 0.25 ML O-passivated iron substrates, are shown. In
both cases, the same dissociative paths determined from equilibrium
calculations at the open surface is displayed. In particular, both
organophosphorus groups separate from the ZDDP molecule, leaving the Zn-S
central units completely isolated. These insights gained by dissociative
reactions are fascinating since no direct experimental evidence exists on the
most favourable ZDDP fragmentation paths. For instance, organophosphorus units
always need to be saturated due to the 5 valence electrons present in the P
atom electronic configuration. Therefore, in the gas phase, these fragments
attracted the S atoms, breaking the Zn-S bonds. On the contrary, the reactive
substrate helps stabilize the P-based groups that leave the molecule and
adsorb on the surface.

Figure 8:  Uniaxial compression optimizations for ZDDP on clean Fe (panel a),
0.25 ML O passivated Fe (panel b) and 1 ML O-passivated Fe (panel c). From
left to right, different snapshots while reaching the target 1 GPa pressure
are reported.

Concerning the O-passivated iron surface at 1 ML, there are no dissociative
reactions even when the external pressure is applied, and the ZDDP molecule
only appears to be stretched. Therefore, the passivation of the iron substrate
substantially reduces the rate of molecular dissociation and at high
O-coverages it impedes it.

##  IV Conclusions

We employed DFT simulation to evaluate key processes involved in the early
stages of ZDDP tribofilm formation, namely: molecular adsorption and its
stress-induced fragmentation both in vacuum and on ferrous substrate with
different level of oxidation.

We found that oxygen passivation significantly reduces the interaction of the
ZDDP molecule with iron. In particular, the energy gain for adsorption
decreases by about one order of magnitude when the Fe (110) surface is fully
oxidised. Moreover, the ZDDP dissociation is energetically favoured over the
clean iron surface and this process is impeded by full oxidation of the iron
substrate.

We studied both uniaxial compression and stretching of ZDDP and we observed
similar fragmentation paths for both cases, forming P-terminated and Zn-S
groups, which serve as the building blocks for the zinc sulfide and phosphate
layers detected in post-mortem analysis [1, 13]. Our findings shed light on
the initial stages of tribofilm formation highlighting the crucial role of the
substrate in favouring molecular dissociation. ZDDP dissociation in oil is
expected to be highly endothermic. However, it has sometimes been proposed as
the key mechanism.

Overall, our results provide a first atomistic insight into ZDDP dissociation,
which represents the first stage of film formation. In particular, while shear
stress application was commonly believed to induce the formation of ZDDP
tribofilms, we showed how the chemical interaction between the molecule and
the surface and how it is influenced, for instance by oxygen, might indeed
play a crucial role as well.

###### Acknowledgements.

These results are part of the SLIDE project, which received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program. (Grant Agreement No. 865633). The authors
want to thank Dr. Margherita Marsili for fruitful discussions.

## References

  * Spikes [2004] H. Spikes, The history and mechanisms of ZDDP, Tribology Letters 17, 469 (2004).
  * Barnes  _et al._ [2001] A. M. Barnes, K. D. Bartle, and V. R. Thibon, A review of zinc dialkyldithiophosphates (zddps): characterisation and role in the lubricating oil, Tribology International 34, 389 (2001).
  * Nicholls  _et al._ [2005] M. A. Nicholls, T. Do, P. R. Norton, M. Kasrai, and G. M. Bancroft, Review of the lubrication of metallic surfaces by zinc dialkyl-dithiophosphates, Tribology international 38, 15 (2005).
  * Gosvami  _et al._ [2015] N. Gosvami, J. Bares, F. Mangolini, A. Konicek, D. Yablon, and R. Carpick, Mechanisms of antiwear tribofilm growth revealed in situ by single-asperity sliding contacts, Science 348, 102 (2015).
  * Neville  _et al._ [2007] A. Neville, A. Morina, T. Haque, and M. Voong, Compatibility between tribological surfaces and lubricant additives—how friction and wear reduction can be controlled by surface/lube synergies, Tribology international 40, 1680 (2007).
  * Martin [1999] J. M. Martin, Antiwear mechanisms of zinc dithiophosphate: a chemical hardness approach, Tribology letters 6, 1 (1999).
  * Zhang and Spikes [2016] J. Zhang and H. Spikes, On the mechanism of ZDDP antiwear film formation, Tribology Letters 63, 24 (2016).
  * Zhang  _et al._ [2020] J. Zhang, J. P. Ewen, M. Ueda, J. S. Wong, and H. A. Spikes, Mechanochemistry of zinc dialkyldithiophosphate on steel surfaces under elastohydrodynamic lubrication conditions, ACS applied materials & interfaces 12, 6662 (2020).
  * Parsaeian  _et al._ [2017] P. Parsaeian, A. Ghanbarzadeh, M. C. Van Eijk, I. Nedelcu, A. Neville, and A. Morina, A new insight into the interfacial mechanisms of the tribofilm formed by zinc dialkyl dithiophosphate, Applied Surface Science 403, 472 (2017).
  * Topolovec-Miklozic  _et al._ [2007] K. Topolovec-Miklozic, T. R. Forbus, and H. A. Spikes, Film thickness and roughness of ZDDP antiwear films, Tribology Letters 26, 161 (2007).
  * Ueda  _et al._ [2019] M. Ueda, A. Kadiric, and H. Spikes, On the crystallinity and durability of ZDDP tribofilm, Tribology Letters 67, 1 (2019).
  * Fang  _et al._ [2023] L. Fang, S. Korres, W. A. Lamberti, M. Webster, and R. W. Carpick, What stress components drive mechanochemistry? A study of ZDDP tribofilm formation, Faraday Discussions 241, 394 (2023).
  * Dorgham  _et al._ [2018] A. Dorgham, A. Azam, A. Morina, and A. Neville, On the transient decomposition and reaction kinetics of zinc dialkyldithiophosphate, ACS Applied Materials & Interfaces 10, 44803 (2018).
  * Heuberger  _et al._ [2007] R. Heuberger, A. Rossi, and N. D. Spencer, XPS study of the influence of temperature on ZnDTP tribofilm composition, Tribology Letters 25, 185 (2007).
  * Sukhomlinov and Müser [2021] S. V. Sukhomlinov and M. H. Müser, Stress anisotropy severely affects zinc phosphate network formation, Tribology Letters 69, 89 (2021).
  * Kumar  _et al._ [2020] D. Kumar, J. Jain, and N. N. Gosvami, In situ study of role of microstructure on antiwear tribofilm formation on az91 magnesium alloy under zinc dialkyldithiophosphate containing lubricant, Advanced Engineering Materials 22, 2000335 (2020), https://doi.org/10.1002/adem.202000335.
  * Mosey  _et al._ [2005] N. J. Mosey, M. H. Muser, and T. K. Woo, Molecular mechanisms for the functionality of lubricant additives, Science 307, 1612 (2005).
  * Chen  _et al._ [2022] Y. Chen, P. Renner, and H. Liang, A review of current understanding in tribochemical reactions involving lubricant additives, Friction 11, 1 (2022).
  * Peeters  _et al._ [2020] S. Peeters, P. Restuccia, S. Loehlé, B. Thiebaut, and M. C. Righi, Tribochemical reactions of MoDTC lubricant additives with iron by Quantum Mechanics/Molecular Mechanics simulations, The Journal of Physical Chemistry C 124, 13688 (2020), doi: 10.1021/acs.jpcc.0c02211.
  * Spikes [2018] H. Spikes, Stress-augmented thermal activation: Tribology feels the force, Friction 6, 1 (2018).
  * Suominen Fuller  _et al._ [2000] M. Suominen Fuller, K. De Jong, M. Kasrai, and G. Bancroft, Electroless generation of phosphate films on metals from zinc dialkyldithiophosphates, Chemistry of materials 12, 1300 (2000).
  * Mee and Torrance [1988] M. Mee and A. Torrance, A study of the lubricated wear of plain bearing materials in a zddp-containing lubricant, Wear 128, 201 (1988).
  * Sheasby  _et al._ [1996] J. Sheasby, T. Caughlin, and W. Mackwood, A comparison of the boundary lubrication of 52100 steel, zirconia and silicon nitride by s, p, s/p and zinc dialkyl dithiophosphate additives, Wear 196, 100 (1996).
  * Mingwu  _et al._ [1993] B. Mingwu, Z. Xushou, and Q. Shangkui, Tribological properties of silicon nitride ceramics coated with molybdenum films under boundary lubrication, Wear 169, 181 (1993).
  * Haque  _et al._ [2007] T. Haque, A. Morina, A. Neville, R. Kapadia, and S. Arrowsmith, Non-ferrous coating/lubricant interactions in tribological contacts: assessment of tribofilms, Tribology International 40, 1603 (2007).
  * Haque  _et al._ [2008] T. Haque, A. Morina, A. Neville, R. Kapadia, and S. Arrawsmith, Study of the zddp antiwear tribofilm formed on the dlc coating using afm and xps techniques, in  _Automotive lubricant testing and advanced additive development_ (ASTM International, 2008).
  * Losi  _et al._ [2021] G. Losi, S. Peeters, F. Delayens, H. Vezin, S. Loehlé, B. Thiebaut, and M. C. Righi, Experimental and ab initio characterization of mononuclear molybdenum dithiocarbamates in lubricant mixtures, Langmuir 37, 4836 (2021).
  * Peeters  _et al._ [2021] S. Peeters, C. Charrin, I. Duron, S. Loehlé, B. Thiebaut, and M. C. Righi, Importance of the catalytic effect of the substrate in the functionality of lubricant additives: the case of molybdenum dithiocarbamates, Materials Today Chemistry 21, 100487 (2021).
  * Ayestarán Latorre  _et al._ [2021] C. Ayestarán Latorre, J. E. Remias, J. D. Moore, H. A. Spikes, D. Dini, and J. P. Ewen, Mechanochemistry of phosphate esters confined between sliding iron surfaces, Communications Chemistry 4, 178 (2021).
  * Mosey and Woo [2003] N. J. Mosey and T. K. Woo, Finite temperature structure and dynamics of zinc dialkyldithiophosphate wear inhibitors: A density functional theory and ab initio molecular dynamics study, The Journal of Physical Chemistry A 107, 5058 (2003), doi: 10.1021/jp034085c.
  * Sukhomlinov  _et al._ [2022] S. V. Sukhomlinov, G. Kickelbick, and M. H. Müser, Mechanochemical ionization: Differentiating pressure-, shear-, and temperature-induced reactions in a model phosphate, Tribology Letters 70, 102 (2022).
  * Mosey and Woo [2004] N. J. Mosey and T. K. Woo, A quantum chemical study of the unimolecular decomposition mechanisms of zinc dialkyldithiophosphate antiwear additives, The Journal of Physical Chemistry A 108, 6001 (2004).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, _et al._ , Quantum espresso: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, _et al._ , Advanced capabilities for materials modelling with quantum espresso, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Giannozzi  _et al._ [2020] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. D. Gironcoli, P. Delugas, F. F. Ruffino, _et al._ , Quantum espresso toward the exascale, The Journal of Chemical Physics 152, 154105 (2020).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Physical Review Letters 77, 3865 (1996).
  * Rappe  _et al._ [1990] A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos, Optimized pseudopotentials, Phys. Rev. B 41, 1227 (1990).
  * Peeters  _et al._ [2022] S. Peeters, A. Barlini, J. Jain, N. N. Gosvami, and M. C. Righi, Adsorption and decomposition of ZDDP on lightweight metallic substrates: Ab initio and experimental insights, Applied Surface Science 600, 153947 (2022).
  * Ta  _et al._ [2023] H. T. T. Ta, N. V. Tran, and M. C. Righi, Atomistic wear mechanisms in diamond: Effects of surface orientation, stress, and interaction with adsorbed molecules (2023), submitted.
  * Pedretti  _et al._ [2023] E. Pedretti, P. Restuccia, and M. C. Righi, Xsorb: a software for identifying the most stable adsorption configuration and energy of a molecule on a crystal surface, Computer Physics Communications , 108827 (2023).
  * Dacre and Bovington [1982] B. Dacre and C. Bovington, The adsorption and desorption of zinc di-isopropyldithiophosphate on steel, Asle Transactions 25, 546 (1982).
  * Tuszynski  _et al._ [2002] W. Tuszynski, J. Molenda, and M. Makowska, Tribochemical conversions of zinc dialkyldithiophosphate (ZDDP) under extremely different pressure conditions, Tribology Letters 13, 103 (2002).
  * Righi  _et al._ [2015] M. C. Righi, S. Loehlé, M. I. de Barros Bouchet, D. Philippon, and J. M. Martin, Trimethyl-phosphite dissociative adsorption on iron by combined first-principle calculations and xps experiments, RSC Advances 5, 101162 (2015).