# A dynamical evolution of GRB-Afterglows in a new generic model

E Zouaoui and N Mebarki  Laboratoire de Physique Mathématique et Subatomique,
Frères Mentouri University, Constantine 1, Algeria<EMAIL_ADDRESS>
###### Abstract

A new model using a general expression of the radiation energy and explaining
the dynamics of the afterglows is proposed. It is shown that this model is
suitable for the ultra-relativistic and non-relativistic phases as well as the
study of radiative and adiabatic fireballs.

##  1 Introduction

Untill now, the mystery of the gamma ray bursts (GRBs) which are the brightest
and violent events in the universe after the Bing Bang is still a big puzzle
for astrophysicists and observers to be understood [1]. In fact, This
phenomenon starts by a prompt emission of the GRBs where the energy $E\sim
10^{51}-10^{54}ergs$, than end up with remnant emission called afterglows. At
this stage scientists were able to measure the redshift [2] [3] from the
obtained spectral lines, showing their cosmological origin. So far, many works
were done in the literature trying to simulate this phenomenon and understand
it [4] [5] [6] [7] [8] [9] [10] [11]. To be more specific the duality of
fireball blast-wave model was successful where its predictions converge with
data. In this model two jets explode from the core in two opposite directions,
as a result they shock the surrounding medium and emit radiations where the
wave lengths run from the gamma to the radio rays spectrum. In what follows we
propose a GRB-afterglows hydrodynamical model taking into account the
synchrotron emission as a major radiation mechanism and neglecting the
absorption (such as synchrotron self-absorption) and diffusion (such as
inverse Compton scattering) effects. It is very important to mention that the
limitation of the earlier models [25] is related to the compatibility with the
Sedov solution starting from a non-relativistic phase where we have an
adiabatic processes and surrounding interstellar medium (ISM). In the present
paper, as it will be shown in the next sections, the Sedov solution was
satisfied achieving the energy conservation in the adiabatic regime. In
section 2, we present the conventional models [12] [13] [15] then we introduce
our generic new model. In section 3, we discuss our numerical results and
finally in section 4, we draw our conclusions.

##  2 Dynamics and radiation

###  2.1 Conventional dynamical models

So far, in literature many models have been proposed to describe the expansion
of the GRB remnants such as the one of Chiang et al of ref.[12] where:

$\frac{d\Gamma}{dm}=-\frac{\Gamma^{2}-1}{M}$ (1)

and

$M=M_{0}+m+(1-\varepsilon)\frac{U}{c^{2}}$ (2)

here $\Gamma$, M, m, $\varepsilon$ and $dU=(\Gamma-1)dmc^{2}$ are the Lorentz
factor, the total mass of the fireball that includes the initial mass $M_{0}$
of the jet (ejected by the progenitor of the GRB), the swept mass from the
external environment, the efficiency of radiation and the fraction of the
internal non radiative energy, respectively. We remind that this first
proposed model was written just for the relativistic phase expansion, because
the solution of eq.(1) is not consistent with Sedov results for an adiabatic
expansion of non relativistic phase. On the other hand, Huang et al of
ref.[13] have proposed new generalized model which is more or less fairly good
for ultra relativistic and non relativistic phases. In this ref.[13], the
authors have used an infinitesimal difference equation for the internal non
radiative energy such as $dU=d[(\Gamma-1)mc^{2}]$ which was assumed and
proposed by Panaitescu et al [14]. In this model, eq.(1) is replaced by the
following equation:

$\frac{d\Gamma}{dm}=-\frac{\Gamma^{2}-1}{M_{0}+\varepsilon
m+2(1-\varepsilon)\Gamma m}$ (3)

where the radiation efficiency $\varepsilon$ was assumed to be constant during
the deceleration.

In 2002 Feng et al [15] [16] suggested that $\varepsilon$ is a variable
belonging to the interval [0,1] and proposed an infinitesimal difference
equation for U of the form:

$dU=d[(1-\varepsilon)U_{ex}]$ (4)

where $U_{ex}=(\Gamma-1)mc^{2}$ is the internal energy produced in the
expansion. In this case one has [17]:

$\frac{d\Gamma}{dm}=-\frac{\Gamma^{2}-1}{M_{0}+m+U/c^{2}+(1-\varepsilon)\Gamma
m}$ (5)

with:

$\varepsilon=\varepsilon_{e}\frac{{t^{{}^{\prime}}}^{-1}_{syn}}{{t^{{}^{\prime}}}^{-1}_{syn}+{t^{{}^{\prime}}}^{-1}_{ex}}$
(6)

Here

$t^{{}^{\prime}}_{syn}=6\pi m_{e}c/(\sigma_{T}B^{{}^{\prime}2}\gamma_{min})$
(7)

and

$t^{{}^{\prime}}_{ex}=R/(\Gamma c)$ (8)

are the synchrotron cooling and expansion times in the co-moving frame
respectively. The parameters $\varepsilon_{e}$ is a fraction of the internal
energy carried by the accelerated electrons in the jet, $m_{e}$ the rest
electron mass, $\sigma_{T}$ is the Thomson cross section and
$B^{{}^{\prime}}=(8\pi\varepsilon_{B}e{{}^{\prime}}){{}^{1/2}}$ the magnetic
energy density where $\varepsilon_{B}$ is the fraction of the kinetic energy
of the shock converted into a magnetic energy and $e{{}^{\prime}}$ an energy
density defined in ref.[18]. The radius R of the blast-wave is determined by
the relation [12]:

$\frac{dR}{dt}=\beta c\Gamma(\Gamma+\sqrt{\Gamma^{2}-1})$ (9)

where $\beta=v/c$ is the jet velocity and $\gamma_{min}$ the minimum Lorentz
factor depending on the index p ($2<p<3$) and is given by [18]:

$\gamma_{min}=\varepsilon_{e}\bigg{(}\frac{p-2}{p-1}\bigg{)}\bigg{(}\frac{m_{p}}{m_{e}}\bigg{)}(\Gamma-1)+1$
(10)

$m_{p}$ being the rest proton mass.

In the highly radiative case where $(\varepsilon\simeq 1,U=0)$, eqs.(3) and
(5), together with eq.(1) reduce to the following differential equation [13]:

$\frac{d\Gamma}{dm}=-\frac{\Gamma^{2}-1}{M_{0}+m}$ (11)

where the analytic solution is [13]:

$\frac{(\Gamma-1)(\Gamma_{0}+1)}{(\Gamma+1)(\Gamma_{0}-1)}=\frac{(m_{0}+M_{0})}{(m+M_{0})}$
(12)

with $\Gamma_{0}$ and $m_{0}$ are the initial values of $\Gamma$ and $m$
respectively.

In the fully adiabatic case where $(\varepsilon\simeq 0,U=U_{ex})$, eqs.(3)
and (5) become [13]:

$\frac{d\Gamma}{dm}=-\frac{\Gamma^{2}-1}{M_{0}+2\Gamma m}$ (13)

leading to the analytic solution [13]:

$(\Gamma-1)M_{0}c^{2}+(\Gamma^{2}-1)mc^{2}=E_{k0}$ (14)

where$E_{k0}$ is the initial value of the kinetic energy $E_{k}$.

###  2.2 Our generic new model

In the Feng et al model [15], and because of the behavior of the radiation
efficiency $\varepsilon$ evolving with time, a differential equation for the
internal residual energy in the fireball for the dynamical evolution of the
GRB-aftergolws was proposed. However the change in the radiative part was
ignored.

In fact the total kinetic energy $E_{k}$ of the fireball is given by [14]:

$E_{k}=(\Gamma-1)(M_{0}+m)+\Gamma U$ (15)

The global conservation of energy implies that:

$dE_{k}=-dE_{rad}$ (16)

where $dE_{rad}$ is the radiative part from the internal energy of the
fireball. The definition given by Blandford et al is
$dE_{rad}=\varepsilon\Gamma(\Gamma-1)dmc^{2}$ [19]. The internal energy
created within the fireball is $U_{ex}$, such that:

$dE_{rad}=\varepsilon\Gamma dU_{ex}$ (17)

This due to the fact that the internal energy can be radiative or remains in
the fireball and has the same form up to a factor $\varepsilon$ (resp.
$1-\varepsilon$) for radiative (resp. residual) internal energy. Therefore,
eq.(17) can be rewritten as:

$dE_{rad}=\varepsilon\Gamma[(\Gamma-1)dm+md\Gamma]c^{2}$ (18)

Using eqs. [(15)- (18)], one obtains the following differential equation
describing the evolution of the fireball and consistent with the Sedov
solution.

$\frac{d\Gamma}{dm}=-\frac{\Gamma^{2}-1}{M_{0}+m(\Gamma+1)+U/c^{2}}$ (19)

It is important to notice that in the adiabatic case when $\varepsilon\simeq
0$, eq.(19) reduces to eq.(13) while in the highly radiative case where
$(\varepsilon\simeq 1)$ eq. (19) gets the new form:

$\frac{d\Gamma}{dm}=-\frac{\Gamma^{2}-1}{M_{0}+m(1+\Gamma)}$ (20)

leading to the following analytic solution:

$\Gamma=-1+\frac{M_{0}}{m}productLog\bigg{[}\frac{m}{M_{0}}e^{\frac{2m-C}{M_{0}}}\bigg{]}$
(21)

(C is an integration constant)

Moreover, it was show that the total luminosity of the fireball has the
following form:

$L=L_{B}+\Gamma\varepsilon mc^{2}\frac{d\Gamma}{dt}$ (22)

where $L_{B}$ is the total luminosity used in the models of refs.[15] [13] and
has as expression $L_{B}=\Gamma\varepsilon[(\Gamma-1)c^{2}dm/dt]$ and $\Gamma$
is giving by eq. (21).

###  2.3 Synchrotron radiation

The distribution of the accelerated electrons by the external shock behind the
blast wave in the absence of the radiation loss is generally assumed to be a
power-law function of the electron energy [17]:

$N^{\prime}_{e}(\gamma_{e})=\frac{dN^{\prime}_{e}}{d\gamma_{e}}=C^{te}\gamma^{-P}_{e},\gamma_{min}\leq\gamma_{e}\leq\gamma_{max}$
(23)

where $\gamma_{max}=a10^{7}(B^{{}^{\prime}}/1G)^{-1/2}$ ref.[19] is the
maximum Lorentz factor and $"a"$ is a factor taking its values between 1 and
10. The power for synchrotron radiation $P_{\nu^{\prime}}$ is defined as [18]:

$P_{\nu^{\prime}}=\frac{2\sqrt{3}e^{2}\nu_{L}}{c}\int^{\gamma_{max}}_{\gamma_{min}}N^{\prime}_{e}(\gamma_{e})F\bigg{(}\frac{\nu^{\prime}}{\nu^{\prime}_{c}}\bigg{)}d\gamma_{e}$
(24)

where F(x) is the synchrotron function such that [20]:

$F(x)=x\int^{\infty}_{x}K_{5/3}(x^{{}^{\prime}})dx^{{}^{\prime}}$ (25)

Here $K_{5/3}$ is the second kind modified Bessel function. The synchrotron
frequency given $\nu^{\prime}_{c}$ by [22]:

$\nu^{\prime}_{c}=\frac{2}{3}\nu_{L}\gamma^{2}_{e}$ (26)

where $\nu_{L}$ is the Lamor frequency:

$\nu_{L}=\frac{1}{2\pi}\frac{eB^{\prime}}{2m_{e}c}$ (27)

Notice that the radiation in the lab frame can be calculated using the
relativistic transformations [21] [22]:

$\nu=\frac{(1+\beta)\Gamma}{1+z}\nu^{{}^{\prime}}$ (28)

where z is the redshift.

The instantenous intensity $F_{\nu}$ giving the light curves at a frequency
$\nu$ in Jansky ($erg.s^{-1}.cm^{-3}.Hz^{-1}$) is:

$F_{\nu}=\frac{1}{4\pi D_{L}(z)^{2}}4\pi\frac{dP_{\nu}}{d\Gamma}$ (29)

where $D_{L}(z)$ is the luminosity distance in the $\Lambda CDM$ model with
$\Omega_{M}=0.3$,$\Omega_{\Lambda}=0.7$ and $H_{0}=71km.s^{-1}Mpc^{-1}$. Note
that for observations by satellites or terrestrial telescopes one has to use
$S_{B}$ by integrating $F_{\nu}$ in a giving interval $[\nu_{1},\nu_{2}]$:

$S_{B}=\int^{\nu_{2}}_{\nu_{2}}F_{\nu}d\nu$ (30)

##  3 Numerical results and discussion

The variation of the mass swepted by the fireball during the deceleration is
given by:

$dm=4\pi R^{2}nm_{p}dR$ (31)

where

$dR=\beta c\Gamma(\Gamma+\sqrt{\Gamma^{2}-1})dt$ (32)

here t is the observed time and measured by telescopes and satellite. In
solving eqs. (1), (3), (5) and (19) numerically, the initial parameters are
taken as: $\Gamma_{0}=250,M_{0}=2.10^{-6}M_{\odot},p=2.1,\theta=10^{\circ}$
(solid angle parameters), $a=4.0$ and $g=0$ with $k=0,n=1cm^{-1}$
corresponding to the interstellar medium (ISM) ref.[25].

Fig.1 displays the evolution of the jet velocity calculated according to eqs.
(3), (5) and (19) [23]. Notice that the three models are compatible with the
Sedov solution, in the non relativistic (NR) limit for the adiabatic case
[24], where the velocity is proportional to $R^{-3/2}$. However in the other
regions the solution of eq.(19) is different from those of eqs. (3), (5). In
fact for the ultra relativistic (UR) case where $x=m/M_{0}<<1$, the principal
branch of Lambert productlog(z) function denoted by $W_{0}(z)$ with
$z=xe^{2x-C}$ (C is integration constant) has the following series expansion
near $\frac{1}{e}$ (e: exponential) as
$W_{0}(x)\approx-1+\sqrt{2(ex+1)}+o(\sqrt{x+\frac{1}{e}})$ which leads to
$\Gamma\approx\frac{-1+e}{\sqrt{2}}+\frac{\sqrt{2}-1}{x}+\mathcal{O}(\sqrt{x+\frac{1}{e}})\sim\frac{\sqrt{2}-1}{x}\sim
R^{-3}$. Of course in this limit $\beta\rightarrow 1$ one can explain the
almost horizontal curve in the region of interval $R\in[10^{14}-10^{18}]cm$,
although there is a small shift due to the fact that in our model
$\beta=1-\frac{1}{\Gamma^{2}}\approx 1\frac{2x}{\sqrt{2}-1}$. For the non
relativistic case where $x>>1$ or equivalently $z>>1$ or
$R\in[10^{19}-10^{23}]cm$ in the region of an interest the $W_{0}(z)$ function
has the following series expansion $W_{0}(z)\approx lnz-
lnlnz+\frac{lnlnz}{lnz}+o(\frac{lnlnz}{lnz})$
$\Gamma=-1+\frac{M_{0}}{m}productLog[\frac{m}{M_{0}}e^{2\frac{m}{M_{0}}-C}]$.
Notice that in our model if $x>>$, $\Gamma\approx
1-\frac{c+ln2}{x}\sim-\frac{c+ln2}{x}\sim R^{-3}$ (numerical result shows that
$c<<-ln2$) leading to $\beta\sim R^{-\frac{3}{2}}$ instead of $\beta\sim
R^{-3}$ ref.[13]-[15]. This explains the difference slope ( in logarithmic
scale between our model and that of ref. [13]-[15] In fact in our new model
the slope is smaller than that of ref.[13]-[15]).

Similar behaviors are noticed in fig.2 where the evolution of the Lorentz
factor for the radiative case shows a faster drop (large slope) in our new
model. Furthermore the curve of the three mentioned models coincide in the
non-relativistic region. In fact in fig.2 (see the curve in UR limit), in the
adiabatic case where $\epsilon=0$ and in all models we obtain the same
differential eq.(13) which has as a solution eq.(14) and therefore in the UR
(resp NR) limit $\Gamma\sim R^{-\frac{3}{2}}$ (resp $\beta\sim
R^{-\frac{3}{2}}$) since $m\sim R^{3}$. For the radiative case where
$\epsilon=1$, in our model we end up with the differential equation
$\frac{d\Gamma}{dm}=-\frac{\Gamma^{2}-1}{M_{0}+m(\Gamma+1)}$ which has as a
solution
$\Gamma=-1+\frac{M_{0}}{m}productLog[\frac{m}{M_{0}}e^{2\frac{m}{M_{0}}-C}]$.
In the UR limit
$\Gamma\approx(\sqrt{2}-1)\frac{M_{0}}{m}+\frac{e-1}{\sqrt{2}}\sim\frac{M_{0}}{m}\sim
R^{-3}$. However in the NR limit one gets $\Gamma\sim
1-(C+ln2)\frac{M_{0}}{m}$ ($\sim 1$) and
$\beta\sim\sqrt{-2\frac{C+ln2}{x}}\bigg{/}(1-\frac{C+ln2}{x})\sim
x^{-\frac{1}{2}}\sim R^{-\frac{3}{2}}$.

It is worth to mention that, in the other models of ref. [13]-[11], one has a
different differential equation that is
$\frac{d\Gamma}{dm}=-\frac{(\Gamma^{2}-1)}{M_{0}+m}$ which has as a solution
$\frac{(\Gamma-1)(\Gamma_{0}+1)}{(\Gamma+1)(\Gamma_{0}-1)}=(\frac{m_{0}+M_{0}}{m+M_{0}})^{2}$
where $m_{0}$ is the initial mass. It is obvious that in the UR limit
$\Gamma\sim\frac{M_{0}}{m}-1\sim R^{-3}$ and in the NR limit $\beta\sim
R^{-3}$ (difference in power of R in comparison to our model). Again, the
difference in the slope of $\Gamma$ as a function of R (in a logarithmic scale
) is due essentially to the fact that we do not have the same solution of the
two different differential equations mentioned above. Finally, in the NR
limit, $\Gamma\rightarrow 1$ that is why all curves tends to the same limit.

Fig.3 shows the total kinetic energy $E_{k}$ as a function of R. Notice that
for the case $\varepsilon\simeq 1.$ $E_{k}$ decreases slower in our new model
than the proposed others models [13] [15] and it is almost flat in the range
of $R\simeq 10^{16}-10^{17}cm$ (UR) and $R\simeq 10^{18}-10^{19}cm$ (NR). In
the adiabatic case $E_{k}$ is constant for all the models $(E_{k}\sim
10^{51}erg)$.

Fig.3 displays the evolution of the total kinetic energy $E_{k}$ as a function
of the distance R in a logarithmic scale for both the radiation and adiabatic
cases where $\varepsilon=1$ and $\varepsilon=0$ respectively. In fact we
remaind that when $\varepsilon=0$ for all models (including ours), the
solution of eq.(13) is that given by eq.(14) leading to
$E_{k}=E_{k_{0}}=(\Gamma-1)M_{0}c^{2}+(\Gamma^{2}-1)mc^{2}=C^{te}$ This
explains why the kinetic energy is constant in dependently of the variation of
R (both in the NR or UR cases). However, for $\varepsilon=1$ (radiative case),
the situation is different,especially in the NR case. The reason lies in the
fact that eq. (21) reduces to $\Gamma\sim 1-(C+ln2)\frac{M_{0}}{m}$ (as it was
pointed out earlier) than the expansion of the kinetic energy
$E_{k}=(\Gamma-1)M_{0}c^{2}$ (common to all models when $\varepsilon=1$) will
be $\sim-(C+ln2)\frac{M}{m}(M+M_{0})$ and since $m<<\gamma_{0}M_{0}$ than
$E_{k}\sim-(C+ln2)=C^{te}$. For the other models, using eq.(12) we deduce that
$E_{k}\sim 2\ \frac{(m_{0}+M_{0})^{2}(\Gamma_{0}-1)}{(\Gamma_{0}+1)}$
$\frac{1}{m+M_{0}}\sim\frac{1}{m}\sim R^{-3}$. In the logarithmic scale, this
behaves like a straight line with a slope $=-3$. For the UR case, our model
gives using eq.(12) and
$\Gamma\sim\frac{\sqrt{2}-1}{m}M_{0}\sim\frac{1}{m}\sim R^{-3}$,
$E_{k}\sim\frac{\Gamma}{m+M_{0}}\sim\frac{\Gamma}{M_{0}}\sim R^{-3}$. Similar
expressions can be obtained for the other models. Notice in this case we have
the same behavior of $E_{k}$ as a function of R expect of course the
proportionality factor a small shift due the difference of the solution of
$\Gamma$ between our model and the others.

Figure 1: Evolution of the jet velocity $\beta=v/c$ as a function of the
distance R (in logarithmic scale). For radiative ($\varepsilon=1$) and
adiabatic ( $\varepsilon=0$) cases. Figure 2:  Evolution of the Lorentz factor
$\Gamma$ as a function of the distance R (in logarithmic scale). For radiative
($\varepsilon=1$) and adiabatic ( $\varepsilon=0$) cases. Figure 3:  Evolution
of the total kinetic energy $E_{k}$ as a function of the distance R (in a
logarithmic scale). For radiative ($\varepsilon=1$) and adiabatic
($\varepsilon=0$) cases. Figure 4:  Evolution of the Lorentz factor $\Gamma$
as a function of the distance R (in logarithmic scale). For
$\varepsilon_{e}=1.$ and $\varepsilon_{B}=0.01$. Figure 5:  Evolution of the
radiative efficiency of the fireball $\varepsilon$ as a function of the
distance R (in logarithmic scale). For $\varepsilon_{e}=1.$ and
$\varepsilon_{B}=0.01$. Figure 6:  Evolution of the total luminosity of the
fireball $L_{bf}$ as a function of the distance R (in logarithmic scale). For
$\varepsilon_{e}=1.$ and $\varepsilon_{B}=0.01$. Figure 7:  Evolution of the
total kinetic energy $E_{k}$ as a function of the distance R (in logarithmic
scale). For $\varepsilon_{e}=1.$ and $\varepsilon_{B}=0.01$. Figure 8:  Light
curve of GRB Afterglow (in logarithmic scale) within Huang et al (1999), Feng
et al (2002) and the new models for radio frequency $\nu_{obs}=3.10^{8}Hz$.
Figure 9:  Integrated light curves of the GRB170208B Afterglow observed by the
XRT / Swift telescope and their simulation (in logarithmic scale).

Fig.5 shows the evolution of the radiative efficiency of the fireball
$\varepsilon$ as a function of the distance R for our model (black solid
line), Feng et al model (red dashed line)and Huang et al (blue dashed line)
for $\varepsilon_{e}=1.$ and $\varepsilon_{B}=0.01$. Notice that our curve is
a little bit shifted in the interval $10^{17}-10^{18}cm$, and coincide in the
region where $R>10^{18}cm$. Notice also in comparison with other models the
maximum which appears at $R\sim 5\times 10^{16}cm$. For $R>>$ (NR case) the
three curves coincide (similarly for the UR case). For the intermediate values
of R, the efficiency $\varepsilon$ in our model is bigger than that of the
Feng and Huang models for the same radius R. For example for $R\sim
210^{17}$cm, $\varepsilon_{ourmodel}\approx 0.49$, $\varepsilon_{Huang}\approx
0.37$ and $\varepsilon_{Feng}\approx 0.27$. To justify this behavior, we
notice that using eqs. (7) and (8) defining $t^{{}^{\prime}}_{syn}$ and
$t^{{}^{\prime}}_{ex}$ respectively, eq. (6) takes the form
$\epsilon=\varepsilon_{e}[1D\Gamma
R^{-1}(A\Gamma+B)^{-1}(\Gamma-1)^{-2}]^{-1}$ where A, B and D are constants.
Now in the NR region where R is large and $\Gamma\simeq 1$,
$\varepsilon\rightarrow\varepsilon_{e}=1$ (see our figure). However in the UR
region where $\Gamma$ is larger
$\varepsilon\rightarrow\varepsilon_{e}\frac{AR\Gamma^{2}\Gamma}{AR\Gamma^{2}+D}$
and this is a decreasing function of $R\Gamma^{2}$.Our numerical results show
that $\frac{d\Gamma}{dt}=\frac{d\Gamma}{dR}\frac{dR}{dt}$ is a decreasing
function of R and since $\frac{d\Gamma}{dt}\sim\Gamma^{2}$ (for a larger value
of $\Gamma$ see eq.(9)) then $\frac{d\Gamma}{dR}$ decreases faster that
$\frac{1}{\Gamma^{2}}$ as a function of R and consequently $\varepsilon$ is a
decreasing function of R as it is shown in fig.5. Now, the shift of
$\varepsilon$ in our model compared to the others, lies in the fact that
numerical results show that within the interval of $R\in[10^{16},10^{19}]$,
the Lorentz factor $\Gamma$ of our model is greater than that of ref.
[13]-[15] for a fixed value of R (see figure 4 ) or equivalently for a fixed
value of $\Gamma$. The radius R of our model is greater than that of ref. ref.
[13]-[15]. Therefore if eq.(6) is rewriten as
$\varepsilon=\varepsilon_{e}\frac{1}{D\Gamma/R(A\Gamma+B)(\Gamma-1)^{2}}$,
than it is obvious that if $R_{ourmodel}$ is greater than $R_{othermodels}$
(for a fixed value of $\Gamma$), one has
$\varepsilon_{ourmodel}>\varepsilon_{othermodels}$.

Fig.6 displays the luminosity as a function of the radius R. Notice that there
is a concordance between the three models in the range $10^{15}-10^{18}cm$
than the curve becomes steeper in our model in the range $>10^{19}cm$.

Fig.6 shows the luminosity $L_{bf}$ as a function of R for a variable
efficiency $\varepsilon$. Notice that L is a decreasing function of R and for
relatively smaller values of R (UR case) it is almost of the same order (for
our model and those of ref. [13] [15]) with a small shift difference (ours is
larger). For example for $R=10^{16}cm$, $L_{ourmodel}\approx
1.4410^{49}erg.s^{-1}$ and $L_{Feng}\approx L_{Huang}\approx
910^{48}erg.s^{-1}$. Now, if R increases, $L_{ourmodel}$ decreases faster than
$L_{Feng}$ and $L_{Huang}$ for example for $R=10^{19}cm$, $L_{ourmodel}\approx
10^{30}erg.s^{-1}$ and $L_{Feng}\approx 3.1310^{33}erg.s^{-1}$ and
$L_{Huang}\approx 1.9210^{34}erg.s^{-1}$. Therefore, for larger values of R
(NR case), the discrepancy increases between the results of our model and
those of ref. [13] [15].This is due mainly to the fact that L contains two
compelling terms $L_{1}=\Gamma\varepsilon(\Gamma-1)\frac{dm}{dt}$ which is
positive and $L_{2}=\Gamma\varepsilon m\frac{d\Gamma}{dt}$ (present only in
our model), which is negative. As R increases, $L_{2}$ increases and
compensates $L_{1}$, such that $L=L_{2}+L_{1}$ become smaller in comparison
with the other models where the term $L_{2}$ is absent. Table 1, summarizes
some numerical values illustrating this fact.

| $L_{1}(erg.s^{-1})$ | $L_{2}(erg.s^{-1})$ | $L_{1}+L_{2}(erg.s^{-1})$
---|---|---|---
R(cm) | Huang | Feng | Our model
$1.09\times 10^{16}$ | $3.48785\times 10^{48}$ | $3.48786\times 10^{48}$ | $3.48798\times 10^{49}$ | $-1.5353\times 10^{46}$ | $3.47263\times 10^{48}$
$1.13\times 10^{17}$ | $5.5\times 10^{47}$ | $4.29\times 10^{47}$ | $2.09\times 10^{48}$ | $-1.16\times 10^{48}$ | $9.36\times 10^{47}$
$1.03\times 10^{18}$ | $2.5\times 10^{41}$ | $2.0\times 10^{40}$ | $1.71\times 10^{41}$ | $-1.40\times 10^{41}$ | $3.08\times 10^{40}$
$1.06\times 10^{19}$ | $1.92\times 10^{34}$ | $3.13\times 10^{33}$ | $1.3323\times 10^{34}$ | $-1.3318\times 10^{34}$ | $5.01\times 10^{30}$
$1.10\times 10^{20}$ | $4.29\times 10^{29}$ | $7.63\times 10^{28}$ | $3.03813\times 10^{29}$ | $-3.03812\times 10^{29}$ | $1.04\times 10^{23}$
Table 1: Some numerical results illustrating the luminosity of the three
models

Fig.7 displays the total kinetic energy $E_{k}$ as a function of the distance
R. For $\varepsilon_{e}=1.$ and $\varepsilon_{B}=0.01$. Notice the discrepancy
between the three models for the range $>5\times 10^{16}cm$. The curves in the
three models became almost flat in the region $R>10^{18}cm$. Notice also the
existence of tow slopes in our new model in the region $R>5\times
10^{16}-10^{18}cm$.

Fig.7 displays $E_{k}$ as a function of R. Notice that for values of
$R=10^{16}cm$ (UR region), $E_{k}$ takes almost the same value for the three
models ($E_{k}\approx 8.86\times 10^{50}erg$). The reason lies to the fact
that in the UR region where $\Gamma$ is larger and $M_{0}\gg m$,
$\varepsilon\approx\varepsilon_{e}$ (see eq.(6)) and therefore in all models
we have $E_{k}\sim(\Gamma-1)M_{0}$ and almost the same differential equation
$\frac{d\Gamma}{dm}=-\frac{\Gamma^{2}-1}{M_{0}}$. In the NR region where $m\gg
M_{0}$ and $\Gamma\approx 1$, $E_{k}$ becomes constant. For example for
$R\approx 10^{19}cm$, $E_{k(Huang)}\approx 1.16\times 10^{50}erg$,
$E_{k(Feng)}\approx 6.09\times 10^{49}erg$ and $E_{k(Ourmodel)}\approx
1.07\times 10^{50}erg$ and $R\approx 10^{20}cm$, $E_{k(Huang)}\approx
1.15\times 10^{50}erg$, $E_{k(Feng)}\approx 5.79\times 10^{49}erg$,
$E_{k(Ourmodel)}\approx 1\times 10^{50}erg$, (see table 2). The reason lies to
the fact that, in the NR region where $\Gamma\approx 1$, $\varepsilon\approx
0$ and $U\approx m(\Gamma-1)+c^{te}$, one has
$\frac{d\Gamma}{dm}\sim-\frac{\Gamma^{2}-1}{M_{0}}$ for all models which leads
to $m\sim\frac{1}{\Gamma-1}$ and therefore $E_{k}\approx 2(\Gamma-1)m+C^{te}$
becomes constant. The difference in the constant between the three models is
due to the complicated expression of $\varepsilon\approx\varepsilon_{e}\times
1/[R(A\Gamma+B)(1-\Gamma)^{2}]$ (if $\Gamma\rightarrow 1$, where $\Gamma$ is a
function of m or equivalently R because slightly different our model to an
other. Any way, numerical results shows that in the NR case, $E_{Huang}$
becomes slightly bigger than ours and almost twice that of Feng and al. In the
intermediate region where $R\approx 3\times 10^{17}cm$; $E_{Ourmodel}\approx
1.68\times 10^{50}cm$, $E_{Feng}\approx 8.81\times 10^{49}cm$ and
$E_{Huang}\approx 1.22\times 10^{50}cm$. In fact, the difference between our
model and Feng et al [15] (resp. Huang et al [13]) is due to the difference in
the expression of the used differential equation $\frac{d\Gamma}{dm}$ (see
eqs. (3) and (5)) (resp. to the difference in the experience of $E_{k}$ (see
eq.(15))and U (see eq.(4) and $U=U_{ex}=(\Gamma-1)mc^{2}$). Table 2,
summarizes the behavior of $E_{k}$ in various regions.

|  $E_{k}$ (erg)
---|---
R(cm) | Huang | Feng | Our model
$10^{16}$ | $8.86\times 10^{50}$ | $8.86\times 10^{50}$ | $8.86\times 10^{50}$
$10^{17}$ | $2.63\times 10^{50}$ | $2.59\times 10^{50}$ | $3.40\times 10^{50}$
$3\times 10^{17}$ | $1.22\times 10^{50}$ | $8.81\times 10^{49}$ | $1.68\times 10^{50}$
$10^{18}$ | $1.16\times 10^{50}$ | $5.09\times 10^{49}$ | $1.07\times 10^{50}$
$10^{19}$ | $1.16\times 10^{50}$ | $5.83\times 10^{49}$ | $1.01\times 10^{50}$
$10^{20}$ | $1.15\times 10^{50}$ | $5.79\times 10^{49}$ | $1.00\times 10^{50}$
Table 2: Some numerical results illustrating the kinetic energy of the three models |  $F_{\nu}$ ($\mu Jy$)
---|---
$t_{obs}$(s) | Huang | Feng | Our model
10 | 789.42 | 789.44 | 789.98
43 | 6045 | 6049 | 7333
$1.2\times 10^{3}$ | 1112 | 890 | 3159
$10^{6}$ | 205 | 98 | 207
$10^{8}$ | 6 | 1 | 4
$10^{9}$ | 0.6 | 0.2 | 0.4
$10^{10}$ | 0.05 | 0.01 | 0.04
Table 3: Some numerical results illustrating the instantenos intensity of
light curve of the three models

Fig.8 shows the light curve of GRB Afterglow within all model for a radio
frequency $\nu_{obs}=3.10^{8}Hz$. Notice the same behavior of the three models
with a small shift in comparison with Feng when $R\leq 10^{2}cm$. Notice also
that:

  1. 1)

$F_{\nu}$ is a decreasing function of t as it is expected.

  2. 2)

There are five regions where the slope $\frac{dF_{\nu}}{dt}$ change its
direction. This is due to the complicated behavior of $\frac{dP}{dt}$ as a
function of t.

  3. 3)

There is a maximum (for all models) around $t\sim 43s$ ($F_{\nu(Huang)}\approx
6045erg.s^{-1}.cm^{-2},F_{\nu(Feng)}\approx
6049erg.s^{-1}.cm^{-2},F_{\nu(Ourmodel)}\approx 7333erg.s^{-1}.cm^{-2}$). This
can be explained by the fact that the synchrotron function K(x) has a maximum
around $x\approx 0.3$ and therefore $F_{\nu}$ should has a maximum
corresponding approximately to $t\approx 43s$

  4. 4)

For relatively small values of t (region 1, $t\sim 10s$), one has almost the
same values of $F_{\nu}$ for qll models. That is: $F_{\nu(Huang)}\approx
789.42erg.s^{-1}.cm^{-2},F_{\nu(Feng)}\approx
789.44erg.s^{-1}.cm^{-2},F_{\nu(Ourmodel)}\approx 789.98erg.s^{-1}.cm^{-2}$.

  5. 5)

In the region 2 (where $\frac{dF_{\nu}}{dt}$ change the direction), our
$F_{\nu}$ is above and different from that of Huang et al and Feng et al. For
example for $t\sim 1.2\times 10^{3}s$ $F_{\nu(Huang)}\approx
1112erg.s^{-1}.cm^{-2},F_{\nu(Feng)}\approx
890erg.s^{-1}.cm^{-2},F_{\nu(Ourmodel)}\approx 3159erg.s^{-1}.cm^{-2}$

  6. 6)

In the region 3 where $t\sim 10^{6}s$ and $F_{\nu(Huang)}\approx
250erg.s^{-1}.cm^{-2},F_{\nu(Feng)}\approx
98erg.s^{-1}.cm^{-2}andF_{\nu(Ourmodel)}\approx 207erg.s^{-1}.cm^{-2}$,
$\frac{dF_{\nu}}{dt}$ changes again its direction (decreases).

  7. 7)

In the region 4 where $F_{\nu}$ drops faster (a more pronounced slope is
noticed) $t\sim 10^{8}-10^{9}s$ and
$F_{\nu(Huang)}\in[6,0.6]erg.s^{-1}.cm^{-2},F_{\nu(Feng)}\in[1,0.2]erg.s^{-1}.cm^{-2}andF_{\nu(Ourmodel)}\in[0.4,0.04]erg.s^{-1}.cm^{-2}$.

  8. 8)

In the region 5 where $t\succeq 10^{10}s(R\approx 10^{19}cm)$ and
$F_{\nu(Huang)}\preceq 0.05erg.s^{-1}.cm^{-2},F_{\nu(Feng)}\preceq
0.01erg.s^{-1}.cm^{-2}$ and $F_{\nu(Ourmodel)}\preceq 0.04erg.s^{-1}.cm^{-2}$,
$F_{\nu}$ changes again its slope (increases).

Notice that in all region 1,2 and 3 in our model $F_{\nu}$ is above and almost
the same as that of Huang et al but slightly different (almost a constant)
from that of Feng et al.

Fig.9 displays the light curves of the GRB170208B Afterglow observed by the
XRT / Swift telescope and their simulation. Notice the fairly good agreement
of our model in the X-Ray Afterglows region $t_{obs}\succeq 3\times 10^{2}s$
and this confirm the viability of our model.

##  4 Conclusion

Throughout this paper, a new generic model of the GRB Afterglows which
verifies the Sedove solution and energy conservation in the adiabatic regime
is constructed. We have studied and discussed the evolution of the jet
velocity, Lorentz factor, kinetic energy, efficiency radiation coefficient,
luminosity and instantenous intensity of light curve as a function of the
radius of the blast wave R and$/$or time t. We have considered both the
radiative and adiabatic cases as well as constant and variable efficiency
radiation coefficient $\varepsilon$. The results f our model are compared with
those of Huang et al and Feng et al models. It turns out that the behaviors
and$/$or numerical results are almost the same in the ultra relativistic
region and different globally in the other region (relativistic and non
relativistic). Finally, fairly good concordance with the observational data of
XRT$/$Swift telescope concerning the integrated light curves of the GRB170208B
was obtained in the interval of our interest $t_{obs}\succeq 3\times 10^{2}s$.

##  5 Acknowledgments

We are very grateful to the Algerian Minister of Higher Educations and
Scientific Research and DGRET for the financial support.

## References

## References

  * [1] Klebesadel R W, Strong I B and Olson R A 1973 _Astrophys. J. Lett._ L85 182
  * [2] Costa E, Frontera F, Heise J, Feroci M, in’t Zand J, Fiore F, Cinti M N, … Piro L and Butler R C 1997 _Nature_ 387 783
  * [3] Metzger, M R, Djorgoski S G, Steidel C C, Kulkarni S R and Frail D A 1997 _Nature_ 387 878
  * [4] Shaviv N J and Dar A 1995 _Mon. Not. Roy. Astron. Soc_ 277 287
  * [5] Sari R and Piran T 1997 _Astrophys. J._ 485 270
  * [6] Dermer C D and Mitman K E 1999 _Astrophys. J._ 513 L5
  * [7] Kumar P 1999 _Astrophys. J._ 523 L113
  * [8] Rees M J and Mészáros P 1994 _Astrophys. J. Lett._ 430 L93
  * [9] Blinnikov S I, Novikov I D, Perevodchikova TV and Polnarev A G 1984 _Soviet. Astron. Lett._ 10 177
  * [10] Katz J I 1994 _Astrophys. J._ 432 107
  * [11] Rees J M and MészÁros P 1992 _Mon. Not. Roy. Astron. Soc_ P258 41
  * [12] Chiang J and Dermer C D 1999 _Astrophys. J._ 512 699
  * [13] Huang Y F, Dai Z G, and Lu T 1999 _MNRAS_ 309 513
  * [14] Panaitescu A, MészÁros P and Rees J M 1998 _Astrophys. J._ 503 315
  * [15] Feng J B, Huang Y F, Dai Z G, and Lu T 2002 _Chin. J. Astron. Astrophys_ 2 525
  * [16] Zouaoui E, Fouka m and Ouichaoui S 2015 _Sciences and Technologie_ A 41 71
  * [17] Dai Z G, Huang Y F and Lu T 1999 _Asrophys. J._ 520 634
  * [18] Huang Y F, Dai Z G, and Lu T 1998 _Astronomy Astrophysics_ 336 L68
  * [19] Sari R, Piran T and Narayan R 1998 _Astrophys. J. Lett._ 497 L17
  * [20] Blandford R D and MeKee C F 1976 _Phys. Fluids._ 19 1130
  * [21] Abramowiotz M and Stegum I A 1965 Handbook of Mathematical Function (New York: Dover)
  * [22] Rybicki G B and Lightman A P 1979 Radiative Procsses in Astrophysics (New York: Wiley)
  * [23] Lind K R and Blandford R D 1985 _Astrophys. J._ 295 358
  * [24] Zouaoui E Fouka M and Ouichaoui S 2011 2010 _Proc. Int, Conf on Algerian Astronomy and Astrophysics (Constantine Algeria)_ vol 1444 (American Institute of Physics) p 359
  * [25] Sedov L 1969 Similarity and Dimensional Methods in Mechanics (New York: Academic)


#  Probabilistic Cauchy Functional Equations

Ehsan Azmoodeh111University of Liverpool, Department of Mathematical Sciences,
L69 7ZL Liverpool, United Kingdom; E-mail address:
<EMAIL_ADDRESS>, Noah Beelders222E-mail address:
<EMAIL_ADDRESS>and Yuliya Mishura333Department of Probability,
Statistics and Actuarial Mathematics, Taras Shevchenko National University of
Kyiv, Mälardalen University, E-mail address<EMAIL_ADDRESS>
###### Abstract

In this short note, we introduce probabilistic Cauchy functional equations,
specifically, functional equations of the following form:

$f(X_{1}+X_{2})\stackrel{{\scriptstyle\textnormal{d}}}{{=}}f(X_{1})+f(X_{2}),$

where $X_{1}$ and $X_{2}$ represent two independent identically distributed
real-valued random variables governed by a distribution $\mu$ having
appropriate support on the real line. The symbol $\stackrel{{\scriptstyle
d}}{{=}}$ denotes equality in distribution. When $\mu$ follows an exponential
distribution, we provide sufficient (regularity) conditions on the function
$f$ to ensure that the unique measurable solution to the above equation is
solely linear. Furthermore, we present some partial results in the general
case, establishing a connection to integrated Cauchy functional equations.

Keywords: Cauchy functional equation; Integrated Cauchy functional equation;
Exponential distribution
MSC 2020: 39B22; 60E05

##  1 Introduction and motivation

Functional equations serve as a cornerstone in various mathematical
disciplines, by playing a pivotal role in both pure and applied mathematics.
Among these, the (additive) Cauchy functional equation

$f(x+y)=f(x)+f(y),\quad x,y\in\mathbb{R},$ (1.1)

stands out as a fundamental and intriguing problem that has garnered
significant attention since its inception in the 19th century. In fact, Cauchy
in 1818 proved that any continuous solution of equation (1.1) is linear; that
is, $f(x)=cx$ for some constant $c$. The argument is straightforward: first
verify linearity for the case of rationals and then pass to the limit using
continuity. In fact, the assumption of continuity can be much further relaxed.
Lebesgue [12] showed that any measurable solution to equation (1.1) is linear.
There is a myriad of interesting literature examining varied conditions on $f$
such that the only solution to equation (1.1) is of a linear form, e.g., [5,
7, 8, 9, 15]. The curious reader can find a comprehensive treatment of the
regularity properties concerning functional equations in Járai’s book [11].
Another fundamental aspect of the Cauchy functional equation (1.1) was
discovered by Hamel in 1905 [10]. He showed that without any conditions on
$f$, the axiom of choice could be used to demonstrate that there exist also
nonlinear solutions of (1.1) for which he also found all such solutions. The
functional equation (1) has been generalized or modified in many other
directions. A typical scenario is to change the domain of validity of the
equation; for instance, to assume that $f$ satisfies equation (1.1) only for
pairs $(x,y)$ belonging to a subset of $\mathbb{R}^{2}$. Additional
information regarding equation (1.1), along with numerous supplementary
references, can be explored in works such as those by Aczél [1] and Aczél-
Dhombres [2]. Moreover, the applications of equation (1.1) in diverse
scientific domains are also extensively covered.

In this paper, we are interested in a probabilistic version of equation (1.1).
For the ease of presentation, we only consider the one dimensional case. Let
$\mu$ be an absolutely continuous probability measure (w.r.t. the Lebesgue
measure) on the real line. The instances of discrete distributions are
excluded since the Cauchy additive functional equation on discrete spaces such
as $\mathbb{Z}_{+}$ is solved trivially. Throughout our presentation, we shall
denote by $S(\mu)$ the support of $\mu$ and additionally assume that
$S(\mu)=\mathbb{R},\,\mathbb{R}_{+}$ or $\mathbb{R}_{-}$ otherwise it is
clearly stated. Also, the probability space $(\Omega,\mathcal{F},{\rm
I\kern-1.60004ptP})$ is fixed and all the random objects are assumed to be
defined on it.

######  Definition 1.1 (Probabilistic Cauchy functional equation).

Fix an integer $n\geq 2$. Let $\mu$ be an absolutely continuous probability
measure on the real line. We say that a function $f:S(\mu)\to\mathbb{R}$
satisfies the $n$-summands probabilistic Cauchy functional equation with
respect to $\mu$ whenever

$f(X_{1}+\dots+X_{n})\stackrel{{\scriptstyle\textnormal{d}}}{{=}}f(X_{1})+\dots+f(X_{n}),$
(1.2)

where $X_{1},\dots,X_{n}$ are $n$ independent identically distributed real-
valued random variables whose laws equal to $\mu$. We simply say that $f$
satisfies the probabilistic Cauchy equation with respect to $\mu$ whenever
equation (1.2) holds with $n=2$.

In contrast to the deterministic Cauchy functional equation, adding more
summands in the probabilistic Cauchy functional equation changes its
distributional properties, and hence we let

$\displaystyle\mathbfcal{PCFE}(\mu;n)$
$\displaystyle:=\Big{\\{}f:S(\mu)\to\mathbb{R}\text{ satisfying equation
}\eqref{eq:PCE}\Big{\\}},\text{ and }$ (1.3)
$\displaystyle\mathbfcal{PCFE}(\mu;2)$ $\displaystyle=\mathbfcal{PCFE}(\mu).$

Here, we have to remove a few ambiguities. First, any solution $f$ to equation
(1.2) remains a valid solution under any suitable modification on every null
set (w.r.t Lebesgue measure). Hence, any solution $f$ to equation (1.2)
induces an equivalent solution class $[f]$ (under obvious equivalence relation
$f\sim g$ iff $f(x)=g(x)$ a.e. Lebesgue measure).Hence, the continuity
assumption would be a natural condition that could be imposed on the solution
$f$ to combine all equivalent solutions into a single solution. Second, when
the support $S(\mu)$ is not the entire real line, any solution
$f\in\mathbfcal{PCFE}(\mu;n)$ can have its domain extended to the whole real
line in many different ways with no harm. In such cases, we are solely
interested in the structure of the solutions of the probabilistic Cauchy
equation on the support of the underlying probability measure, and hence
specify the domain of solution $f\in\mathbfcal{PCFE}(\mu;n)$ as $S(\mu)$. A
typical instance is when $\mu\sim\textbf{Exp}(\lambda=1)$ is an exponential
distribution whose support is $\mathbb{R}_{+}$ (we will examine it in more
details later on). Now, we ask ourselves the following natural question
sharing the same spirit of the classical Cauchy function equation (1.1):

What regularity conditions dictate that any solution
$f\in\mathbfcal{PCFE}(\mu;n)$ takes the linear form $f(x)=cx$ (for some
constant $c$) almost everywhere or everywhere?

In this note, we consider the above question with a particular emphasis on
when $\mu$ is an exponential distribution. We address in Section 1.2, one of
the main reasons for choosing the exponential distribution is that the above
central question arose naturally during the study of a useful property in the
class of (fractional) phase-type distributions. In Section 2.1 we present
several sufficient regularity conditions on $f$ that guarantee the linear form
of the unique solution to the probabilistic Cauchy equation (1.2). Section 2.2
compiles a few partial results for a general probability measure $\mu$ with
connections to integrated Cauchy functional equations.

###  1.1 Almost sure Cauchy functional equation

There exist other stochastic versions of the Cauchy functional equation (1.1)
proposed in the existing literature. For example, in [13], the authors studied
solutions of the stochastic Cauchy functional equation

$f(X_{1}+X_{2})=f(X_{1})+f(X_{2}),\quad a.s.,$ (1.4)

where $X_{1},X_{2}\sim\mathcal{N}(0,1)$ are two independent standard Gaussian
random variables. Note that, clearly, equation (1.4) is much stronger than our
probabilistic Cauchy functional equation (1.2) when $\mu=\mathcal{N}(0,1)$.
However, their primary goal differs from ours in that they connect the almost
sure Cauchy functional equation (1.4) to a martingale problem. In fact, they
proved the following interesting observation. Let $W=(W_{t})$ be a Brownian
motion w.r.t. its natural filtration. Then, the stochastic process
$M_{t}=f(W_{t})$ is a martingale where $f$ is a solution to equation (1.4) iff
$f(x)=cx$ almost everywhere w.r.t. Lebesgue measure. Additionally, if one
imposes the extra regularity assumption of ”right-continuity” on the
martingale $M$, then the ”almost everywhere” condition in the former statement
can be removed. Their study was inspired by the work of Smirnov [16] who
showed that a weak version of Bernstein’s characterization of the Gaussian
distribution implies the local integrability of a measurable solution of the
Cauchy functional equation (1.1).

###  1.2 Connection with the inhomogeneous phase-type distributions

Let’s consider the parameterised exponential family of distributions; that is

$\mathbfcal{E}:=\\{F(\theta;x):=1-e^{-\theta
x}\,:\,\theta>0,x\in(0,\infty)\\}.$ (1.5)

It is a classical result (or can be seen by direct computation) that the
exponential family $\mathbfcal{E}$ enjoys the following interesting property:
for $n\in\mathbb{N}$ and every $\theta\in(0,\infty)$,

$\underbrace{F(\theta;\cdot)\star\cdots\star F(\theta;\cdot)}_{n\text{ times
}}=F^{\star
n}(\theta;\cdot)=\frac{(-\theta)^{n}}{(n-1)!}\,\frac{\partial^{n-1}}{\partial\theta^{n-1}}G(\theta;\cdot),$
(1.6)

where $G(\theta;\cdot):=-\theta^{-1}F(\theta;\cdot)$ and $\star$ stands for
the convolution operator. The beauty of relation (1.6) is that it reduces the
derivation of the distribution of the sum of $n$ independent exponential
random variables to simple differentiation w.r.t. the parameter variable. Now,
let $X(\theta)\sim F(\theta;\cdot)$ be a random variable corresponding to the
distribution function $F(\theta;\cdot)$ for each $\theta\in\Theta$. In
addition, let $h:\mathbb{R}\rightarrow\mathbb{R}$ be a bijective function, and
consider the new random variable

$Y(\theta):=h(X(\theta))\sim
F_{h}(\theta;\cdot):=F\big{(}\theta;h^{-1}(\cdot)\big{)},\quad\forall\,\theta\in\Theta.$
(1.7)

Denote by
$\mathbfcal{E}_{h}:=\\{F_{h}(\theta;\cdot)=F(\theta;h^{-1}(\cdot))\,:\,\theta\in\Theta\\}$
the parameterised family of probability distributions corresponding to random
variables $Y(\theta)$, $\theta\in\Theta$. Then, in light of relation (1.7),
family $\mathbfscr{E}_{h}$ naturally inherits property (1.6) from the
exponential family $\mathbfscr{E}$. More precisely, for every natural number
$n$, we define the $h$-convolution operator

$Y(\theta)^{\mathbin{\mathchoice{\ooalign{$\displaystyle\vbox{\hbox{\scalebox{0.77778}{$\displaystyle\bigcirc$}}}$\cr$\displaystyle\star$\cr}}{\ooalign{$\textstyle\vbox{\hbox{\scalebox{0.77778}{$\textstyle\bigcirc$}}}$\cr$\textstyle\star$\cr}}{\ooalign{$\scriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptstyle\bigcirc$}}}$\cr$\scriptstyle\star$\cr}}{\ooalign{$\scriptscriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptscriptstyle\bigcirc$}}}$\cr$\scriptscriptstyle\star$\cr}}}n}:=h\left(X_{1}(\theta)+\dots+X_{n}(\theta)\right),\quad\theta\in\Theta,$
(1.8)

where the random variables $X_{i}(\theta)\sim F(\theta;\cdot)$ are
independent. We also adapt the conventional notation
$F_{h}(\theta;\cdot)^{\mathbin{\mathchoice{\ooalign{$\displaystyle\vbox{\hbox{\scalebox{0.77778}{$\displaystyle\bigcirc$}}}$\cr$\displaystyle\star$\cr}}{\ooalign{$\textstyle\vbox{\hbox{\scalebox{0.77778}{$\textstyle\bigcirc$}}}$\cr$\textstyle\star$\cr}}{\ooalign{$\scriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptstyle\bigcirc$}}}$\cr$\scriptstyle\star$\cr}}{\ooalign{$\scriptscriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptscriptstyle\bigcirc$}}}$\cr$\scriptscriptstyle\star$\cr}}}n}$
for the distribution of the random variable
$Y(\theta)^{\mathbin{\mathchoice{\ooalign{$\displaystyle\vbox{\hbox{\scalebox{0.77778}{$\displaystyle\bigcirc$}}}$\cr$\displaystyle\star$\cr}}{\ooalign{$\textstyle\vbox{\hbox{\scalebox{0.77778}{$\textstyle\bigcirc$}}}$\cr$\textstyle\star$\cr}}{\ooalign{$\scriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptstyle\bigcirc$}}}$\cr$\scriptstyle\star$\cr}}{\ooalign{$\scriptscriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptscriptstyle\bigcirc$}}}$\cr$\scriptscriptstyle\star$\cr}}}n}$.
Now, it is straightforward to use property (1.6) to find distributional
properties of
$Y(\theta)^{\mathbin{\mathchoice{\ooalign{$\displaystyle\vbox{\hbox{\scalebox{0.77778}{$\displaystyle\bigcirc$}}}$\cr$\displaystyle\star$\cr}}{\ooalign{$\textstyle\vbox{\hbox{\scalebox{0.77778}{$\textstyle\bigcirc$}}}$\cr$\textstyle\star$\cr}}{\ooalign{$\scriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptstyle\bigcirc$}}}$\cr$\scriptstyle\star$\cr}}{\ooalign{$\scriptscriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptscriptstyle\bigcirc$}}}$\cr$\scriptscriptstyle\star$\cr}}}n}$.
In this vein, note that, for $n\geq 2$,

$\displaystyle
F_{h}(\theta;x)^{\mathbin{\mathchoice{\ooalign{$\displaystyle\vbox{\hbox{\scalebox{0.77778}{$\displaystyle\bigcirc$}}}$\cr$\displaystyle\star$\cr}}{\ooalign{$\textstyle\vbox{\hbox{\scalebox{0.77778}{$\textstyle\bigcirc$}}}$\cr$\textstyle\star$\cr}}{\ooalign{$\scriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptstyle\bigcirc$}}}$\cr$\scriptstyle\star$\cr}}{\ooalign{$\scriptscriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptscriptstyle\bigcirc$}}}$\cr$\scriptscriptstyle\star$\cr}}}n}$
$\displaystyle:={\rm
I\kern-1.60004ptP}\left(Y(\theta)^{\mathbin{\mathchoice{\ooalign{$\displaystyle\vbox{\hbox{\scalebox{0.77778}{$\displaystyle\bigcirc$}}}$\cr$\displaystyle\star$\cr}}{\ooalign{$\textstyle\vbox{\hbox{\scalebox{0.77778}{$\textstyle\bigcirc$}}}$\cr$\textstyle\star$\cr}}{\ooalign{$\scriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptstyle\bigcirc$}}}$\cr$\scriptstyle\star$\cr}}{\ooalign{$\scriptscriptstyle\vbox{\hbox{\scalebox{0.77778}{$\scriptscriptstyle\bigcirc$}}}$\cr$\scriptscriptstyle\star$\cr}}}n}\leq
x\right)={\rm
I\kern-1.60004ptP}\big{(}h(X_{1}(\theta)+\dots+X_{n}(\theta))\leq
x\big{)}=F\big{(}\theta;h^{-1}(x)\big{)}^{\star n}$ (1.25)
$\displaystyle=\frac{(-\theta)^{n}}{(n-1)!}\frac{\partial^{n-1}}{\partial\theta^{n-1}}G(\theta;h^{-1}(x)).$
(1.26)

It is mentioned that the distributions of the form
$Y=h(X_{1}(\theta)+\dots+X_{n}(\theta))$ bears an intimate relation to the
inhomogeneous phase-type matrix distributions presented in [3]. Hence, it is
meaningful, due to the validity of property (1.25), to ask which functions $h$
permit that

$h\left(X_{1}(\theta)+\dots+X_{n}(\theta)\right)\stackrel{{\scriptstyle\textnormal{d}}}{{=}}h(X_{1}(\theta))+\dots+h(X_{n}(\theta)).$
(1.27)

##  2 Main results

###  2.1 $\mathbfcal{PCFE}(\mu)$: case of the exponential distribution

######  Theorem 2.1.

Let the following assumptions be fulfilled.

  * $(i)$

The independent random variables $X_{1},X_{2}\sim\mu$, where $\mu$ is an
exponential distribution with parameter $\lambda>0$.

  * $(ii)$

The real-valued function $f$ has the following properties: $f\in
C^{1}[0,+\infty),f(0)=0$, $f(+\infty)=+\infty$, $f^{\prime}(0)>0$ and $f$ is
strictly increasing.

  * $(iii)$

$f\in\mathbfcal{PCFE}(\mu)$, i.e., the following random variables are
identically distributed:

$f\left(X_{1}+X_{2}\right)\stackrel{{\scriptstyle\scriptsize\textnormal{d}}}{{=}}f\left(X_{1}\right)+f\left(X_{2}\right).$
(2.1)

Then $f(x)=cx$ for some $c\geqslant 0$.

######  Remark 2.1.

  * $(a)$

Assuming that $f\in C^{1}[0,+\infty)$, we consider the right-hand derivative
at zero, as usual.

  * $(b)$

If we assume statements $(i)$ and $(iii)$ from the above theorem as well as
that $f$ is strictly increasing and $f\in C^{1}[0,+\infty)$, then
automatically $f(0)=0$ and $f(+\infty)=+\infty$ otherwise (2.1) is obviously
violated. Indeed, if, for example, $f(+\infty)=f_{\infty}\in\mathbb{R}$, then
${\rm
I\kern-1.60004ptP}\left(f\left(X_{1}\right)+f\left(X_{2}\right)>f_{\infty}\right)>0$
whilst ${\rm
I\kern-1.60004ptP}\left(f\left(X_{1}+X_{2}\right)>f_{\infty}\right)=0$.

  * $(c)$

We can consider symmetric assumption: $f\in C^{1}[0,+\infty),f(0)=0$,
$f(+\infty)=-\infty$, and $f$ is strictly decreasing. Then we shall get that
$f(x)=cx$ for some $c\leqslant 0$, for which the proof will be the same, but
with obvious symmetric changes.

###### Proof.

Note that (2.1) is equivalent to condition ${\rm
I\kern-1.60004ptP}(A_{x})={\rm I\kern-1.60004ptP}(B_{x})$ for any
$x\in\mathbb{R}$, where

$A_{x}=\left\\{f\left(X_{1}+X_{2}\right)\geqslant x\right\\},\quad
B_{x}=\left\\{f\left(X_{1}\right)+f\left(X_{2}\right)\geqslant
x\right\\}\text{. }$

Consider these events separately. Since $f(0)=0,f(+\infty)=+\infty$ and $f$ is
strictly increasing, there exists the inverse function $f^{-1}$ that is also
strictly increasing on $\mathbb{R}^{+}=[0,+\infty)$, $f^{-1}(0)=0$,
$f^{-1}(+\infty)=+\infty$. Then for any $x>0$

$\displaystyle A_{x}=$ $\displaystyle\;\left\\{X_{1}+X_{2}\geqslant
f^{-1}(x)\right\\}$ $\displaystyle=$ $\displaystyle\;\left\\{X_{2}\geqslant
f^{-1}(x)\right\\}\cup\left\\{0<X_{2}<f^{-1}(x),X_{1}\geqslant
f^{-1}(x)-X_{2}\right\\}\text{, }$

and

$\displaystyle B_{x}=$ $\displaystyle\;\left\\{f\left(X_{1}\right)\geqslant
x-f\left(X_{2}\right)\right\\}$ $\displaystyle=$
$\displaystyle\;\left\\{X_{2}\geqslant
f^{-1}(x)\right\\}\cup\left\\{X_{2}<f^{-1}(x),X_{1}\geqslant
f^{-1}\left(x-f\left(X_{2}\right)\right)\right\\}.$

Therefore,

${\rm I\kern-1.60004ptP}\left(A_{x}\right)={\rm
I\kern-1.60004ptP}\left(B_{x}\right)$ if and only if

${\rm I\kern-1.60004ptP}\left\\{0<X_{2}<f^{-1}(x),X_{1}\geqslant
f^{-1}(x)-X_{2}\right\\}={\rm
I\kern-1.60004ptP}\\{0<X_{2}<f^{-1}(x),X_{1}\geqslant
f^{-1}(x-f\left(X_{2}\right))\\}.$ (2.2)

In terms of exponential densities, Eq. (2.2) can be rewritten in the following
equivalent form

$\int_{0}^{f^{-1}(x)}\lambda e^{-\lambda u}\int_{f^{-1}(x)-u}^{\infty}\lambda
e^{-\lambda v}\textnormal{d}v\textnormal{d}u=\int_{0}^{f^{-1}(x)}\lambda
e^{-\lambda u}\int_{f^{-1}(x-f(u))}^{\infty}\lambda e^{-\lambda
v}\textnormal{d}v\textnormal{d}u,$ (2.3)

for any $x>0$. In turn, it follows from equation (2.3) that

$\int_{0}^{f^{-1}(x)}\textnormal{exp}\\{-\lambda
u-\lambda(f^{-1}(x)-u)\\}\textnormal{d}u=\int_{0}^{f^{-1}(x)}\textnormal{exp}\\{-\lambda
u-\lambda f^{-1}(x-f(u))\\}\textnormal{d}u,$

and hence that

$f^{-1}(x)e^{-\lambda
f^{-1}(x)}=\int_{0}^{f^{-1}(x)}\textnormal{exp}\\{-\lambda u-\lambda
f^{-1}(x-f(u))\\}\textnormal{d}u,$ (2.4)

for any $x>0$, because

$-\lambda u-\lambda(f^{-1}(x)-u)=-\lambda f^{-1}(x)$

and does not depend on $u$.

Now, since $f^{-1}(x)$ can be an arbitrary number $y\geq 0$, we can rewrite
Eq. (2.4) as

$e^{-\lambda y}\cdot y=\int_{0}^{y}\exp\\{-\lambda u-\lambda
f^{-1}(f(y)-f(u))\\}\textnormal{d}u,$

for any $y\geq 0$. Furthermore, by implementing the change of variables
$y-u=v$ under the sign of the right-hand side integral, the latter can be
rewritten as

$e^{-\lambda y}\cdot y=\int_{0}^{y}\exp\left\\{-\lambda(y-v)-\lambda
f^{-1}(f(y)-f(y-v))\right\\}\textnormal{d}v,$

or, equivalently as

$y=\int_{0}^{y}\exp\left\\{\lambda v-\lambda
f^{-1}(f(y)-f(y-v))\right\\}\textnormal{d}v.$ (2.5)

As a consequence of the equation (2.5), we can state that

$\psi(y):=\int_{0}^{y}\left(\exp\left\\{\lambda
v-f^{-1}(f(y)-f(y-v))\right\\}-1\right)\textnormal{d}v=0,\quad y\geq 0,$ (2.6)

which will be our function of interest.

Since we assume that $f\in C^{1}[0,+\infty)$ and since $\psi(0)=0$, Eq. (2.6)
is equivalent to the following equation: $\psi^{\prime}(y)=0$. Now, we take
explicitly the derivative of the right-hand side of Eq. (2.6) and take into
account that

$\exp\left\\{\lambda y-\lambda
f^{-1}(f(y)-f(y-y))\right\\}-1=\exp\left\\{\lambda y-\lambda
f^{-1}(f(y))\right\\}-1=0.$

Therefore, by recalling for continuously differentiable functions that

$\partial_{y}\Bigl{(}\int_{0}^{y}\varphi(y,v)\textnormal{d}v\Bigr{)}=\int_{0}^{y}\partial_{y}\varphi(y,v)\textnormal{d}v+\varphi(y,y),$

and recalling also the formula for the differentiation of the inverse
function, we get that

$\psi^{\prime}(y)=\int_{0}^{y}\exp\\{\theta(y,v)\\}\times\frac{f^{\prime}(y)-f^{\prime}(y-v)}{f^{\prime}(f^{-1}(f(y)-f(y-v)))}\textnormal{d}v=0,$
(2.7)

where

$\theta(y,v)=\lambda v-\lambda f^{-1}(f(y)-f(y-v)).$

Note that $f^{\prime}(x)>0$ for any $x>0$ because the function $f$ is strictly
increasing, and $f^{\prime}(x)>0$ by assumption of the theorem. Therefore

$f^{\prime}\left(f^{-1}(f(y)-f(y-v))\right)>0,$

for any $v\in[0,y]$, and all components under the sign of integral are
continuous implying that the integral is correctly defined.

Our goal now is to prove that the numerator $f^{\prime}(y)-f^{\prime}(y-v)$ is
a zero function. Thus, consider the behavior of $f^{\prime}(y)$, for $y>0$. Of
course, there can only be the two following possibilities.

  1. $(i)$

There exists a strictly increasing sequence of points $y_{n}^{up},\,n\geq 1$
such that $y_{n}^{up}\uparrow\infty$ as $n\rightarrow\infty$, and

$f^{\prime}\left(y_{n}^{up}\right)\geq f^{\prime}(y)\quad\text{for all}\quad
y\leq y_{n}^{up},\quad\text{and}\quad
f^{\prime}\left(y_{n+1}^{up}\right)>f^{\prime}\left(y_{n}^{up}\right),\,n\geq
1.$

Since by Eq. (2.7)

$\exp\\{\theta(y,v)\\}>0\quad\text{and}\quad
f^{\prime}(f^{-1}(f(y)-f(y-v)))>0,$

we conclude that $\psi^{\prime}\left(y_{n}^{up}\right)$ can be zero only if
$f^{\prime}(y_{n}^{up})=f^{\prime}(y_{n}^{up}-v)$ for all
$v\in[0,y_{n}^{up}]$; i.e. $f^{\prime}(y)=c_{n}$ for all $y\leq y_{n}^{up}$
and some $c_{n}>0$. Changing $n$ for $n+1$, we conclude that $c_{n}=c_{n+1},$
and so $f^{\prime}(y)=c_{n}=c_{1}$ for all $y\leq y_{n+1}^{up}$.

Finally, since $y_{n}^{up}\uparrow\infty$ as $n\rightarrow\infty$, we get that
$c_{n}=c$ for some $c>0$, and $f^{\prime}(y)=c$ for all $y\geq 0$. Therefore,
$f(y)=cy$ for $y\geq 0$.

Note that the case of symmetric downwards conditions (see Remark 2.1, $(c)$)
with $y_{n}^{down}\uparrow\infty$, where
$f^{\prime}\left(y_{n}^{down}\right)\leq f^{\prime}(y)$ for all $y\leq
y_{n}^{down}$ and
$f^{\prime}\left(y_{n+1}^{down}\right)<f^{\prime}\left(y_{n}^{down}\right)$,
can be considered similarly.

  2. $(ii)$

There exists two points $y_{0}^{up}$ and $y^{down}_{0}$ such that
$f^{\prime}\left(y_{0}^{up}\right)\geq f^{\prime}(y)$ for all $y\geq 0$ and
$f^{\prime}({y}_{0}^{down})\leq f^{\prime}(y)$ for all $y\geq 0$.

Let $y_{0}^{up}<{y}_{0}^{down}$ (the case ${y}_{0}^{up}={y}_{0}^{{down}}$ is
trivial, and the case ${y}_{0}^{up}>{y}_{0}^{down}$ can be considered
similarly). Then we can ”forget” about $y_{0}^{up}$, consider the whole
interval $[0,{y}_{0}^{down}]$ and deduce exactly as in $(i)$, case $(b)$, that
$f^{\prime}(y)=c$ for all $0\leq y\leq{y}_{0}^{down}$. However, it means
immediately that

$f^{\prime}({y}_{0}^{up})=f^{\prime}({y}_{0}^{down}),$

and therefore $f^{\prime}(y)=c$ for all $y\geq 0.$

∎

######  Remark 2.2.

If it is assumed that $f\in C^{2}[0,\infty)$, then the assumption
$f^{\prime}(0)>0$ can be avoided. Indeed, suppose for some $\alpha\in(0,1)$
that $\lim_{x\rightarrow
0}\frac{f^{\prime}(x)}{x^{\alpha}}=c_{\alpha}\in(0,\infty)$. Then by
l’Hospital’s rule

$\displaystyle\lim_{x\rightarrow
y}\frac{f^{\prime}(y)-f^{\prime}(x)}{f^{\prime}\left(f^{-1}(f(y)-f(x))\right)}=$
$\displaystyle\;c_{\alpha}^{-1}\lim_{x\rightarrow
y}\frac{f^{\prime}(y)-f^{\prime}(x)}{\bigl{[}f^{-1}(f(y)-f(x))\bigr{]}^{\alpha}}$
$\displaystyle=$ $\displaystyle\;(\alpha c_{\alpha})^{-1}\lim_{x\rightarrow
y}\frac{-f^{\prime\prime}(x)\bigl{[}f^{-1}(f(y)-f(x))\bigr{]}^{1-\alpha}}{-\left[f^{\prime}\left(f^{-1}(f(y)-f(x))\right)\right]^{-1}f^{\prime}(x)}\quad$
$\displaystyle=$ $\displaystyle\;(\alpha
c_{\alpha})^{-1}\frac{f^{\prime\prime}(y)}{f^{\prime}(y)}\lim_{x\rightarrow
y}\bigl{[}f^{-1}(f(y)-f(x))\bigr{]}^{1-\alpha}f^{\prime}(f^{-1}(f(y)-f(x)))$
$\displaystyle=$ $\displaystyle\;0,$

for any $y>0$, and where integral (2.7) is correctly defined.

We also observe that the assumption of independence of
$X_{1},X_{2}\sim\textnormal{Exp}(\lambda)$ for $\lambda>0$ yields for any
suitable (e.g. bounded) $g(x_{1},x_{2}):\mathbb{R}^{2}\rightarrow\mathbb{R}$
that

${\rm
I\kern-1.60004ptE}g(X_{1},X_{2})=\int_{\mathbb{R}^{2}}g(x_{1},x_{2})\lambda^{2}\textnormal{e}^{-\lambda(x_{1}+x_{2})}\textnormal{d}x_{1}\textnormal{d}x_{2}.$

Using this, we have the following.

######  Proposition 2.1.

Let $f:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}\in C[0,\infty)$ s.t.
$f(x+y)\geq f(x)+f(y)$ for all $x,y\geq 0$. Assume further that
$f\in\mathbfcal{PCFE}(\mu)$ where $\mu$ is an exponential distribution with
parameter $\lambda>0$. Then $f(x)=cx$ for some $c>0$.

###### Proof.

Let $X_{1},X_{2}\sim\textnormal{Exp}(\lambda)$ be two independent random
variables. Consider

${\rm
I\kern-1.60004ptE}\textnormal{e}^{-f(X_{1}+X_{2})}=\int_{\mathbb{R}^{2}}\textnormal{e}^{-f(x+y)}\lambda^{2}\textnormal{e}^{-\lambda(x+y)}\textnormal{d}x\,\textnormal{d}y,$

and

${\rm
I\kern-1.60004ptE}\textnormal{e}^{-f(X_{1})-f(X_{2})}=\int_{\mathbb{R}^{2}}\textnormal{e}^{-f(x)-f(y)}\lambda^{2}\textnormal{e}^{-\lambda(x+y)}\textnormal{d}x\,\textnormal{d}y,$

Then ${\rm I\kern-1.60004ptE}\textnormal{e}^{-f(X_{1}+X_{2})}={\rm
I\kern-1.60004ptE}\textnormal{e}^{-f(X_{1})-f(X_{2})}$ by the assumption that
$f(X_{1}+X_{2})\stackrel{{\scriptstyle\scriptsize\textnormal{d}}}{{=}}f(X_{1})+f(X_{2})$,
but

$\int_{\mathbb{R}^{2}}\textnormal{e}^{-f(x+y)-\lambda x-\lambda
y}\textnormal{d}x\,\textnormal{d}y\leq\int_{\mathbb{R}^{2}}\textnormal{e}^{-f(x)-f(y)-\lambda
x-\lambda y}\textnormal{d}x\,\textnormal{d}y,$

because the $\textnormal{e}^{-f(x+y)-\lambda x-\lambda
y}\leq\textnormal{e}^{-f(x)-f(y)-\lambda x-\lambda y}$ for any
$(x,y)\in\mathbb{R}^{2}$. Both statements can only simultaneously be true if
$f(x+y)=f(x)+f(y)$, for all $x,y$. Consequently, by the continuity of $f$, we
have that $f(x)=cx$ for some $c>0$. ∎

######  Remark 2.3.

The case where $f:\mathbb{R}_{+}\rightarrow(-\infty,0]$ and $f(x+y)\leq
f(x)+f(y)$ can be considered similarly.

We conclude this section by presenting two results regarding general
probability measures, leveraging the assumption of sub(super)additivity within
the solution function $f$.

######  Proposition 2.2.

Let $f\in\mathbfscr{PCFE}(\mu)$ be a continuous solution to the probabilistic
Cauchy functional equation (1.2). Assume further that
$f:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is strictly increasing, $f(0)=0$,
$f(+\infty)=+\infty$ and that $f(x+y)\leq(\geq)f(x)+f(y)$ for all $x,y$. Then
$f(x)=cx$ for some constant $c$.

###### Proof.

Let $p$ denote the density of $\mu$. In this case, we observe that Eq. (2.3)
takes the form

$\int_{0}^{f^{-1}(x)}p(u)\int_{f^{-1}(x)-u}^{\infty}p(v)\textnormal{d}v\textnormal{d}u=\int_{0}^{f^{-1}(x)}p(u)\int_{f^{-1}(x-f(u))}^{\infty}p(v)\textnormal{d}v\textnormal{d}u.$

Denoting $f^{-1}(x)=y\geq 0$, we get from the above equation that

$\int_{0}^{y}p(u)\int_{y-u}^{\infty}p(v)\textnormal{d}v\textnormal{d}u=\int_{0}^{y}p(u)\int_{f^{-1}(f(y)-f(u))}^{\infty}p(v)\textnormal{d}v\textnormal{d}u,$

or equivalently that

$\int_{0}^{y}p(u)\int^{f^{-1}(f(y)-f(u))}_{y-u}p(v)\textnormal{d}v\textnormal{d}u=0.$
(2.8)

Assuming that $f(x+y)\geq f(x)+f(y)$ for all $x,y\geq 0$, we observe that
$f(y-u)\leq f(y)-f(u)$ for all $0\leq u\leq y$, whence $y-u\leq
f^{-1}\bigl{(}f(y)-f(u)\bigr{)}$ for all $0\leq u\leq y$. Therefore, Eq. (2.8)
can hold only for the case when $f(y-u)=f(y)-f(u)$ for all $0\leq u\leq y$. ∎

######  Proposition 2.3.

Let $n\geq 2$ be an integer and $f\in\mathbfscr{PCFE}(\mu;n)$ be a continuous
solution to the probabilistic Cauchy functional equation (1.2) such that ${\rm
I\kern-1.60004ptE}_{\mu}|f|<\infty$. Assume further that
$f(x+y)\geq(\leq)f(x)+f(y)$ for all $x,y$. Then $f(x)=cx$ for some constant
$c$.

###### Proof.

Let $p$ denote the density of $\mu$. Then, the validity of equation (1.2), the
integrability assumption and a direct conditioning argument yield that

$\int_{S(\mu)^{n}}\big{\\{}f(x_{1}+\dots+x_{n})-\left(f(x_{1})+\dots+f(x_{n})\right)\big{\\}}p(x_{1})\dots
p(x_{n})\textnormal{d}x_{1}\dots\textnormal{d}x_{n}=0.$ (2.9)

Now, the assumption of super(sub)-linearity together with (2.9) implies that
the function

$g:S(\mu)^{n}\to\mathbb{R};\,(x_{1},\dots,x_{n})\mapsto
f(x_{1}+\dots+x_{n})-\left(f(x_{1})+\dots+f(x_{n})\right)$ (2.10)

is zero almost everywhere for all $(x_{1},\dots,x_{n})\in
S(\mu)^{n}(=\mathbb{R}^{n},\mathbb{R}^{n}_{+},\text{ or
}\mathbb{R}^{n}_{-}\text{ by our assumption})$. But, note that $g$ is a
jointly continuous function, and, hence, it is identically zero everywhere.
Finally, using the continuity of $f$ along with the simple fact that $f(0)=0$,
we immediately conclude that $f(x+y)=f(x)+f(y)$ for all $(x,y)$, and thus the
claim follows at once by employing the classic Cauchy functional equation. ∎

###  2.2 $\mathbfcal{PCFE}(\mu)$: a partial result in the general case and its
connection with integrated Cauchy functional equations

In this section, we provide a partial result in the general case beyond the
exponential distribution via tools available in the realm of integrated Cauchy
functional equations (ICFEs). In order to briefly motivate the latter, let
$X\sim\mu$ be a random variable distributed as $\mu$ being an absolutely
continuous probability measure on the real line with density $p$. Under the
integrability assumption ${\rm I\kern-1.60004ptE}|f(X)|<\infty$, one can
easily see that the probabilistic Cauchy functional equation (1.2) (equating
the first moment of both sides with $n=2)$ yields that

$\int_{S(\mu)}\int_{S(\mu)}\big{\\{}f(x+y)-\left(f(x)+f(y)\right)\big{\\}}p(x)p(y)\textnormal{d}x\textnormal{d}y=0.$
(2.11)

Note that, obviously, any integrable solution $f\in\mathbfcal{PCFE}(\mu)$ is a
solution to integral equation (2.11), however the inverse implication is far
from being true. For example, consider $\mu=\mathcal{N}(0,1)$, a standard
Gaussian distribution. Then, one can easily check that function $f(x)=x^{3}$
is a solution to equation (2.11) but not a solution to the probabilistic
Cauchy functional equation (1.2).

Now observe that a sufficient requirement for the validity of the integral
equation (2.11) is that, for almost all $x$,

$f(x)+{\rm I\kern-1.60004ptE}f(X)=\int_{S(\mu)}f(x+y)\mu(\textnormal{d}y).$
(2.12)

The latter requirement (2.12) may be viewed as an integrated version of the
additive Cauchy functional equation $f(x+y)=f(x)+f(y)$ that appears as a
variant of the so-called Lau-Rao theorem (see Section 2 of [14]). The origin
of the Lau-Rao theorem can be traced back to the celebrated Choquet-Deny
convolution equation $\mu=\mu\ast\sigma$ (see [6] for details). For the sake
of completeness, below, we present a version taken from [14, Section 4],
slightly adapted to fit our framework.

######  Theorem 2.2 (The additive Lau-Rao Theorem).

Let $\mu$ be a non-arithmetic $\sigma$-finite measure on $\mathbb{R}_{+}$ such
that $\mu(\\{0\\})<1$. Consider the following equation for almost all
$x\in\mathbb{R}_{+}$ (w.r.t Lebesgue measure):

$\int_{\mathbb{R}_{+}}f(x+y)\mu(\textnormal{d}y)=f(x)+{\rm
I\kern-1.60004ptE}f(X),$ (2.13)

where $f:\mathbb{R}_{+}\to\mathbb{R}$ is a locally integrable Borel function.
Assume that $f$ is either an increasing or decreasing function almost
everywhere. Then, every solution to the equation (2.13) takes the form

$f(x)=\begin{cases}a+b\left(1-e^{-\eta x}\right)&\mbox{ if }\,\eta\neq 0\\\
cx+d&\mbox{ if }\,\eta=0\end{cases}$ (2.14)

where $a,b,c,d$ are constants and $\eta$ is such that
$\int_{\mathbb{R}_{+}}e^{-\eta x}\mu(\textnormal{d}x)=1$. Furthermore, assume
that all the assumptions above prevail apart from that instead of
$\mathbb{R}_{+}$, we consider the whole real line $\mathbb{R}$. Then, every
solution to the equation (2.13) takes the form
$\lambda_{1}f_{1}+\lambda_{2}f_{2}$ where the non-negative constants
$\lambda_{1}$ and $\lambda_{2}$ are so that $\lambda_{1}+\lambda_{2}=1$ and,
$f_{1}$ and $f_{2}$ are of the form (2.14) with $\eta$ replaced respectively
with $\eta_{1}$ and $\eta_{2}$ satisfying in equation
$\int_{\mathbb{R}}e^{-\eta_{i}x}\mu(\textnormal{d}x)=1$, $i=1,2$.

######  Remark 2.4.

In our case, we have that $\mu$ is a probability distribution with a support
$S(\mu)=\mathbb{R}_{+},\mathbb{R}_{-}$ or $\mathbb{R}$. For
$S(\mu)=\mathbb{R}_{+},\mathbb{R}_{-}$, we have that
$\int_{S(\mu)}\textnormal{e}^{-\eta x}\mu(\textnormal{d}x)=1$ occurs only for
$\eta=0$. However, for $S(\mu)=\mathbb{R}$, we have that
$L_{\mu}(\eta):=\int_{\mathbb{R}}\textnormal{e}^{-\eta
x}\mu(\textnormal{d}x)=1$ for $\eta=0$ and possibly one other $\eta\neq 0$.

Indeed, let
$\mathcal{D}_{\eta}:=\\{\eta\in\mathbb{R}:L_{\mu}(\eta)<\infty\\}$. First note
that $\mathcal{D}_{\eta}$ may be a singleton; for instance, the Cauchy
distribution has $\mathcal{D}_{\eta}=\\{0\\}$. If, however, some $\eta_{u}>0$
is in $\mathcal{D}_{\eta}$, then since we have for $s\in[0,\eta_{u}]$ that
$\textnormal{e}^{-(\eta_{u}-s)x}\leq\textnormal{e}^{-\eta_{u}x}$ on
$x\in\mathbb{R}_{-}$ and $\textnormal{e}^{-sx}\leq 1$ on $x\in\mathbb{R}_{+}$,
it must be that $L_{\mu}(\eta)$ exists at least on the interval
$[0,\eta_{u}]$. A similar reasoning shows that if $\eta_{d}<0$ is in
$\mathcal{D}_{\eta}$, then so is the interval $[\eta_{d},0]$. Therefore,
$\mathcal{D}_{\eta}$ must itself be some interval (not necessarily including
the endpoints).

Furthermore, since the convexity of $\textnormal{e}^{x}$ yields that
$\textnormal{e}^{\lambda
x_{1}+(1-\lambda)x_{2}}\leq\lambda\textnormal{e}^{x_{1}}+(1-\lambda)\textnormal{e}^{x_{2}}$
for $x_{1},x_{2}\in\mathcal{D}_{\eta}$, $\lambda\in[0,1]$, taking expectations
of this inequality yields that $L_{\mu}(\eta)$ is convex. In turn, since
$L_{\mu}(\eta)=1$ occurs for $\eta=0$, it must be that $L_{\mu}(\eta)=1$
occurs only for one other $\eta\neq 0$ at most. Indeed, assuming existence of
two nonzero points, we get from convexity of $L_{\mu}$ the existence of the
whole interval $[0,\eta_{u})$ where ${\rm
I\kern-1.60004ptE}\textnormal{e}^{-\eta X}=1$ for all $\eta\in[0,\eta_{u})$.
Then we can choose $2\eta\in[0,\eta_{u})$ such that ${\rm
I\kern-1.60004ptE}\textnormal{e}^{-2\eta X}=1$. Therefore,
$\text{Var}(\textnormal{e}^{-\eta X})={\rm
I\kern-1.60004ptE}\textnormal{e}^{-2\eta X}-\bigl{(}{\rm
I\kern-1.60004ptE}\textnormal{e}^{-\eta X}\bigr{)}^{2}=0$ showing that $X=0$
a.s. A simple example of existence of two points where $L_{\mu}(\eta)=1$ is
the distribution $\mu\sim\mathcal{N}(\gamma,\sigma^{2})$ with mean $\gamma\neq
0$, variance $\sigma^{2}$, having
$L_{\mu}(\eta)=\exp\\{-\gamma\eta+\sigma^{2}\eta^{2}/2\\}$.

ICFEs have played a significant role in characterisation problems appearing in
non-contemporary probability theory and mathematical statistics. For details,
the reader can consult textbook [4]. The next proposition offers a partial
result to our central question by employing this approach.

######  Theorem 2.3.

Fix an integer $n\geq 2$. Let $\mu$ be a non-degenerate probability measure
with support $S(\mu)$ having finite second moment. Let
$f:S(\mu)(=\mathbb{R}_{+},\mathbb{R}_{-},\text{or }\mathbb{R})\to\mathbb{R}\in
L^{2}(\mu)$ be a Borel measurable function so that with $X\sim\mu$ it holds
that for almost all $x$

$\hskip 113.81102pt{\rm
I\kern-1.60004ptE}\left(f(x+X)-f(X)\right)\geq(\leq)f(x).\hskip
142.26378pt{\mathbf{(H)}}$

Assume that $f\in\mathbfcal{PCFE}(\mu;n)$ be either an increasing or
decreasing function almost everywhere. Then, $f(x)=cx$ for some constant $c$.

###### Proof of Theorem 2.3.

First, let $X_{1},\dots,X_{n-1}$ be $n$ mutually independent random variables
having distribution $\mu$. Then, using the tower property and by repetitively
using assumption (H), we obtain that

$\displaystyle{\rm I\kern-1.60004ptE}$
$\displaystyle\bigl{(}f(x+X_{1}+\cdots+X_{n-1})\bigr{)}={\rm
I\kern-1.60004ptE}\bigl{(}{\rm
I\kern-1.60004ptE}\bigl{(}f(x+X_{1}+\cdots+X_{n-1})\big{|}X_{1},\dots,X_{n-2}\bigr{)}\bigr{)}$
$\displaystyle\geq(\leq){\rm I\kern-1.60004ptE}\Big{(}{\rm
I\kern-1.60004ptE}\bigl{(}f(x+X_{1}+\cdots+X_{n-2})\big{|}X_{1},\dots,X_{n-3}\bigr{)}+f(X_{n-1})\Big{)}$
$\displaystyle\;\;\;\;\;\vdots$ $\displaystyle\geq(\leq){\rm
I\kern-1.60004ptE}\Big{(}{\rm
I\kern-1.60004ptE}\bigl{(}f(x+X_{1}\bigr{)}+f(X_{2})+\cdots+f(X_{n-1})\Big{)}$
$\displaystyle\geq(\leq){\rm
I\kern-1.60004ptE}\bigl{(}f(x)+f(X_{1})+\cdots+f(X_{n-1})\bigr{)}.$

The latter is equivalently expressed as

$\int_{S(\mu)^{n-1}}\Big{(}f(x+x_{1}+\dots+x_{n-1})-\left(f(x)+f(x_{1})+\dots+f(x_{n-1})\right)\Big{)}\mu(\textnormal{d}x_{1})\dots\mu(\textnormal{d}x_{n-1})\geq(\leq)0.$
(2.15)

Now, assume that $f\in\mathbfcal{PCFE}(\mu;n)$ such that $f\in L^{2}(\mu)$. By
equating the expectations of both sides of equation (1.2), we obtain

$\int_{S(\mu)^{n}}\Big{(}f(x+x_{1}+\dots+x_{n-1})-\left(f(x)+f(x_{1})+\dots+f(x_{n-1})\right)\Big{)}\mu(\textnormal{d}x_{1})\dots\mu(\textnormal{d}x_{n-1})\mu(\textnormal{d}x)=0,$

and so the above statement combined with Eq. (2.15) yields for almost all $x$
that

$f(x)+(n-1){\rm
I\kern-1.60004ptE}f(X)=\int_{S(\mu)^{n-1}}f(x+x_{1}+\dots+x_{n-1})\mu(\textnormal{d}x_{1})\dots\mu(\textnormal{d}x_{n-1}).$
(2.16)

Now, for $n\geq 3$, observe from Eq. (2.16) and assumption (H) that

$\displaystyle f(x)+(n-2){\rm I\kern-1.60004ptE}f(X)$ $\displaystyle={\rm
I\kern-1.60004ptE}\bigl{(}f(x+X_{1}+\cdots+X_{n-2}+X_{n-1})-f(X_{n-1})\bigr{)}$
$\displaystyle={\rm I\kern-1.60004ptE}\Big{(}{\rm
I\kern-1.60004ptE}\bigl{(}f(x+X_{1}+\cdots+X_{n-2}+X_{n-1})-f(X_{n-1})\big{|}X_{1},\dots
X_{n-2}\bigr{)}\Big{)}$ $\displaystyle\geq(\leq){\rm
I\kern-1.60004ptE}\bigl{(}f(x+X_{1}+\cdots+X_{n-2})\bigr{)},$

and so repeating the above procedure yields that $f(x)+{\rm
I\kern-1.60004ptE}f(X)\geq(\leq){\rm I\kern-1.60004ptE}\bigl{(}f(x+X)\bigr{)}$
where $X\sim\mu$. Hence, combining the latter inequalities with assumption
(H), we obtain that for almost all $x$,

$f(x)+{\rm I\kern-1.60004ptE}f(X)=\int_{S(\mu)}f(x+y)\mu(\textnormal{d}y).$
(2.17)

Now, we consider two cases separately.
Case $S(\mu)=\mathbb{R}_{+}$: First we note that equation
$L_{\mu}(\eta):=\int_{\mathbb{R}_{+}}e^{-\eta x}\mu(\textnormal{d}x)=1$ has
only trivial solution $\eta=0$. Now, since it is assumed that $f$ is monotone,
a direct application of the Lau-Rao Theorem 2.2 yields that all the solutions
to equation (2.17) are given by $f(x)=cx+d$. Referring to the probabilistic
Cauchy functional equation, we deduce that $d=0$. The case
$S(\mu)=\mathbb{R}_{-}$ follows similarly.

Case $S(\mu)=\mathbb{R}$: In addition to $\eta=0$, there may possibly be an
additional non-trivial $(\eta\neq 0)$ solution to the equation
$L_{\mu}(\eta)=1$, see Remark 2.4.

Relying on the Lau-Rao Theorem 2.2 since $f$ is monotone by assumption, every
solution $f$ to the equation (2.17) can be written as a convex combination
$f=\lambda_{1}f_{1}+\lambda_{2}f_{2}$ of two solutions of the form given by
(2.14). The rest of the proof is dedicated to showing that the only valid
solutions to the probabilistic Cauchy equation (1.2) are given by linear
functions. Hence, assume that $\eta_{1}\neq 0$ and $\eta_{2}=0$ are the
solutions of $L_{\mu}(\eta)=1$, and consider (for $\lambda_{1}\in[0,1]$ and
$\lambda_{2}=1-\lambda_{1}$)

$f(x)=\lambda_{1}b(1-\textnormal{e}^{-\eta_{1}x})+\lambda_{2}cx,$ (2.18)

for some constants $b$ and $c$ that are chosen to ensure the monotonicity of
$f$; for instance, if $\eta<0$, then $b=1$ and $c=-1$ are permissible
constants. (We exclude the the constants $a,d$ in (2.14) since they must equal
$0$ for $f$ to satisfy in the probabilitistic Cauchy functional equation.)

To avoid the trivial cases, let $\lambda_{1}\neq 0,1$. Let
$\lambda:=\lambda_{2}c/(\lambda_{1}b)$. Since $f\in\mathbfcal{PCFE}(\mu;n)$
and $f\in L^{2}(\mu)$ by our assumption, the second moments of
$f(\sum_{i=1}^{n}X_{i})$ and $\sum_{i=1}^{n}f(X_{i})$ exist and must be equal
where $f$ is given by Eq. (2.18); i.e. we have that

$\displaystyle{\rm
I\kern-1.60004ptE}\Bigl{(}1-\exp(-\eta_{1}\sum_{i=1}^{n}X_{i})+\lambda\sum_{i=1}^{n}X_{i})\Bigr{)}^{2}$
$\displaystyle={\rm
I\kern-1.60004ptE}\Bigl{(}\sum_{i=1}^{n}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}+\lambda\sum_{i=1}^{n}X_{i}\Bigr{)}^{2}.$
(2.19)

Expanding the r.h.s. of the above equation, we get

$\displaystyle{\rm
I\kern-1.60004ptE}\Bigl{(}\sum_{i=1}^{n}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}+\lambda\sum_{i=1}^{n}X_{i}\Bigr{)}^{2}$
$\displaystyle={\rm
I\kern-1.60004ptE}\Bigl{(}\sum^{n}_{i=1}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}\Bigr{)}^{2}+\lambda^{2}{\rm
I\kern-1.60004ptE}\Bigl{(}\sum_{i=1}^{n}X_{i}\Bigr{)}^{2}+2\lambda{\rm
I\kern-1.60004ptE}\Bigl{(}\sum^{n}_{i=1}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}\sum_{i=1}^{n}X_{i}\Bigr{)},$
$\displaystyle={\rm
I\kern-1.60004ptE}\Bigl{(}\sum^{n}_{i=1}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}\Bigr{)}^{2}+\lambda^{2}{\rm
I\kern-1.60004ptE}\Bigl{(}\sum_{i=1}^{n}X_{i}\Bigr{)}^{2}+2\lambda{\rm
I\kern-1.60004ptE}\Bigl{(}\sum_{i=1}^{n}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}X_{i}\Bigr{)}.$
(2.20)

and the l.h.s. yields

$\displaystyle{\rm
I\kern-1.60004ptE}\Bigl{(}1-\exp(-\eta_{1}\sum_{i=1}^{n}X_{i})+\lambda\sum_{i=1}^{n}X_{i})\Bigr{)}^{2}$
$\displaystyle={\rm
I\kern-1.60004ptE}\Bigl{(}1-\exp(-\eta_{1}\sum_{i=1}^{n}X_{i})\Bigr{)}^{2}+\lambda^{2}{\rm
I\kern-1.60004ptE}\Bigl{(}\sum_{i=1}^{n}X_{i}\Bigr{)}^{2}+2\lambda{\rm
I\kern-1.60004ptE}\Bigl{(}\bigl{[}1-\exp(-\eta_{1}\sum_{i=1}^{n}X_{i})\bigr{]}\sum_{i=1}^{n}X_{i}\Bigr{)},$
(2.21)

However, by using that $L_{\mu}(\eta_{1})=1$ and the independence assumption
of $X_{1},...,X_{n}$, the last term of Eq. (2.21) becomes

$\displaystyle{\rm
I\kern-1.60004ptE}\Bigl{(}\bigl{[}1-\exp(-\eta_{1}\sum_{i=1}^{n}X_{i})\bigr{]}\sum_{i=1}^{n}X_{i}\Bigr{)}$
$\displaystyle=\sum_{i=1}^{n}{\rm
I\kern-1.60004ptE}\Bigl{(}X_{i}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}\Bigr{)}$
$\displaystyle\;\;\;\;\;+\sum_{i=1}^{n}{\rm
I\kern-1.60004ptE}\Bigl{(}X_{i}\bigl{[}1-\exp(-\eta_{1}\sum_{j\neq
i}X_{j})\bigr{]}\Bigr{)}$ $\displaystyle={\rm
I\kern-1.60004ptE}\Bigl{(}\sum_{i=1}^{n}X_{i}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}\Bigr{)},$

which is the same as that in Eq. (2.20). Therefore, by using the above
identity and equating Eq. (2.21) and (2.20), we have that Eq. (2.19) reduces
to

${\rm
I\kern-1.60004ptE}\Bigl{(}1-\exp(-\eta_{1}\sum^{n}_{i=1}X_{i})\Bigr{)}^{2}={\rm
I\kern-1.60004ptE}\Bigl{(}\sum^{n}_{i=1}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}\Bigr{)}^{2}.$
(2.22)

Since $L_{\mu}(\eta_{1})=1$, we observe for the l.h.s. of the above equation
that

$\displaystyle{\rm
I\kern-1.60004ptE}\Bigl{(}1-\exp(-\eta_{1}\sum_{i=1}^{n}X_{i})\Bigr{)}^{2}$
$\displaystyle={\rm
I\kern-1.60004ptE}\Bigl{(}\exp(-2\eta_{1}X_{1})\Bigr{)}^{n}-1,$ (2.23)

and, by using independence of the $X_{i}$ to notice that the expectations of
the cross-terms are zero, that the r.h.s. yields

$\displaystyle{\rm
I\kern-1.60004ptE}\Bigl{(}\sum^{n}_{i=1}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}\Bigr{)}^{2}$
$\displaystyle=\sum^{n}_{i=1}{\rm
I\kern-1.60004ptE}\Bigl{(}\bigl{[}1-\exp(-\eta_{1}X_{i})\bigr{]}\Bigr{)}^{2}$
$\displaystyle=n\Bigl{(}{\rm
I\kern-1.60004ptE}\Bigl{(}\exp(-2\eta_{1}X_{1})\Bigr{)}-1\Bigr{)}.$ (2.24)

Now, we shall denote $x:={\rm I\kern-1.60004ptE}\exp(-2\eta_{1}X_{1})$, and
thus we have by equating Eqs. (2.23) and (2.24) that

$(x^{n}-1)=n(x-1),$

which is equivalent to

$(x-1)(\sum^{n-1}_{k=0}x^{k}-n)=0.$ (2.25)

If $x<1$, we have that $\sum^{n-1}_{k=0}x^{k}<n$, and if $x>1$, then
$\sum^{n-1}_{k=0}x^{k}>n$, and so
$(x-1)(\sum^{n-1}_{k=0}x^{k}-n)>0$ if $x\neq 1$. Thus Eq. (2.25) yields that
$x=1$. On the other hand, note that $x\neq 1$ yields a contradiction.
Therefore, functions $f$ of the form (2.18) cannot be a solution to the
probabilistic Cauchy equation (1.2).

∎

######  Remark 2.5.

In the proof of Theorem 2.3, when the support $S(\mu)=\mathbb{R}$ is the
entire real line, following the Lau-Rao Theorem 2.2, we had to treat solutions
of the non-linear forms given by (2.18) involving a non-trivial $(\eta\neq 0)$
constant $\eta$ so that $L_{\mu}(\eta)=1$. However, all those non-linear
solutions would be immediately eliminated as soon as measure $\mu$ is
symmetric about the origin, i.e., $\mu(A)=\mu(-A)$ or equivalently
$F(-x)=1-F(x)$ where $F$ stands for the corresponding distribution function.
In other words, the only solution to the equation $L_{\mu}(\eta)=1$ is given
by the trivial solution $\eta=0$. To see this, Theorem 3c from pp. 239 of [17]
guarantees that $L_{\mu}(\eta)=\eta\int_{\mathbb{R}}e^{-\eta
x}F(x)\textnormal{d}x$ for $\eta\neq 0$, and provided that the Laplace
transform $L_{\mu}(\eta)$ exists. Hence, using the symmetry assumption, we
have for $\eta>0$ that

$\displaystyle\frac{L_{\mu}(\eta)}{\eta}$
$\displaystyle=\int_{0}^{\infty}e^{-\eta
x}\bigl{(}F(x)-1+1)\textnormal{d}x+\int_{-\infty}^{0}e^{-\eta
x}F(x)\textnormal{d}x$ $\displaystyle=-\int_{0}^{\infty}e^{-\eta
x}\bigl{(}1-F(x)\bigr{)}\textnormal{d}x+\frac{1}{\eta}+\int^{\infty}_{0}\textnormal{e}^{\eta
x}F(-x)\textnormal{d}x$ $\displaystyle=-\int_{0}^{\infty}e^{-\eta
x}\bigl{(}1-F(x)\bigr{)}\textnormal{d}x+\frac{1}{\eta}+\int^{\infty}_{0}\textnormal{e}^{\eta
x}\bigl{(}1-F(x)\bigr{)}\textnormal{d}x,$

which consequently yields that

$\frac{L_{\mu}(\eta)-1}{\eta}=\int^{\infty}_{0}\bigl{(}\textnormal{e}^{\eta
x}-\textnormal{e}^{-\eta x}\bigr{)}\bigl{(}1-F(x)\bigr{)}\textnormal{d}x,$

a contradiction if $L_{\mu}(\eta)=1$ for $\eta>0$ (the case of $\eta<0$ can be
shown to yield a similar contradiction). A simple example fulfilling this
setup is given by the centred Gaussian distribution
$\mu=\mathcal{N}(0,\sigma^{2})$ with
$L_{\mu}(\eta)=e^{-\sigma^{2}\eta^{2}/2}$.

######  Remark 2.6.

For a given probability measure $\mu$, assumption (H) can be viewed as a form
of super(sub) linearity on average with respect to the measure $\mu$, which,
in general, is much weaker than the classic pointwise super(sub) linearity. In
other words, for a given function $f$ assumption (H) is equivalent to

$\int_{S(\mu)}f(x+y)\mu(\textnormal{d}y)\geq(\leq)\int_{S(\mu)}\left(f(x)+f(y)\right)\mu(\textnormal{d}y),\quad\text{
for almost all $x$}.$ (2.26)

######  Remark 2.7.

This remark aims to clarify the relationships between all the assumptions
imposed on the solution function $f$ in the aforementioned results.

(a) It is clear that any function satisfying in $f(x+y)\geq f(x)+f(y)$
satisfies also assumption (H) no matter the choice of the underlying
probability measure $\mu$. However, the converse is not necessarily true.
Indeed, Lemma 3.1 from the Appendix demonstrates a piecewise linear function
$f$ satisfying assumption (H) w.r.t. the exponential distribution but for
which pointwise super-linearity fails for a set which is a non-trivial
interval.

(b) Let $\mu$ be an exponential distribution. The requirement that
$f:[0,\infty)\mapsto[0,\infty)$ is increasing (decreasing) does not
automatically imply the validity of assumption $\mathbf{(H)}$, and hence both
assumptions are required for Proposition 2.3. For example, consider the
function $f(x)=e^{x/2}$ which is clearly (strictly) increasing but violates
$\mathbf{(H)}$. Furthermore, Lemma 3.1 shows that the converse does not hold
either since a non-increasing function satisfying in assumption $\mathbf{(H)}$
exists by choosing $a=b=c=d=1$ and $r\geq 2e-1$.

(c) Note that for a given function $f:[0,\infty)\to[0,\infty)$ the condition
$f(x+y)\geq f(x)+f(y),\;\forall x,y\in\mathbb{R}_{+}$ implies that $f$ is
necessarily increasing but not strictly increasing (consider $f(x)=0$, for
$x\in[0,1]$ and $f(x)=x-1$ for $x\geq 1$). Furthermore, the strictly
increasing assumption does not imply super-linearity, i.e., $f(x+y)\geq
f(x)+f(y),\;\forall x,y\in\mathbb{R}_{+}$. Indeed, consider $f(x)=2x-\sin x$.
Then, finding $x$ and $y$ s.t. $f(x+y)\leq f(x)+f(y)$ is the same as finding
$x$ and $y$ s.t. $\sin(x+y)\geq\sin x+\sin y$, and so choosing
$x,y\in(\pi,\tfrac{3}{2}\pi)$ is sufficient.

## References

  * [1] Aczél, J. Lectures on Functional Equations and Their Applications.  Academic Press, Jan. 1966.
  * [2] Aczél, J., and Dhombres, J. Functional Equations in Several Variables.  Cambridge University Press, 1989.
  * [3] Albrecher, H., and Bladt, M. Inhomogeneous phase-type distributions and heavy tails.  J. Appl. Probab. 56, 4 (2019), 1044–1064.
  * [4] Balasubrahmanyan, R., and Lau, K.-S. Functional Equations in Probability Theory.  Elsevier, 2014.
  * [5] Banach, S. Sur l’équation fonctionnelle f(x+y)=f(x)+f(y).  Fundamenta Mathematicae 1, 1 (1920), 123–124.
  * [6] Choquet, G., and Deny, J. Sur l’équation de convolution $\mu=\mu\ast\sigma$.  C. R. Acad. Sci. Paris 250 (1960), 799–801.
  * [7] Darboux, G. Sur la composition des forces en statique.  Bulletin des sciences mathématiques et astronomiques 9 (1875), 281–288.
  * [8] Darboux, M. G. Sur le théorème fondamental de la géométrie projective.  Mathematische Annalen 17, 1 (1880), 55–61.
  * [9] Fréchet, M. Pri la funkcia equacio f(x+y)=f(x)+f(y).  Enseignement des Mathematiques (1913).
  * [10] Hamel, G. Eine Basis aller Zahlen und die unstetigen Lösungen der Funktionalgleichung: f(x+y)=f(x)+f(y).  Mathematische Annalen 60, 3 (1905), 459–462.
  * [11] Járai, A. Regularity Properties of Functional Equations in Several Variables.  Springer Science & Business Media, 2005.
  * [12] Lebesgue, H. Sur les transformations ponctuelles.  Atti della R. Accademia delle Scienze di Torino 42 (1907).
  * [13] Mania, M., and Tikanadze, L. Functional equations and martingales.  Aequationes Mathematicae 96, 1 (2022), 221–241.
  * [14] Radhakrishna Rao, C., and Shanbhag, D. N. Recent results on characterization of probability distributions: a unified approach through extensions of Deny’s theorem.  Advances in Applied Probability 18, 3 (1986), 660–678.
  * [15] Sierpiński, W. Sur l’équation fonctionnelle f(x+y)=f(x)+f(y).  Fundamenta Mathematicae 1, 1 (1920), 116–122.
  * [16] Smirnov, S. N. A probabilistic note on the Cauchy functional equation.  Aequationes Mathematicae 93, 2 (2019), 445–449.
  * [17] Widder, D. The Laplace Transform.  Princeton University Press, 1941.

##  3 Appendix

######  Lemma 3.1.

Let $X$ be an exponentially distributed random variable with parameter
$\lambda=1$. Then, there exist constants $a,b,c,d,r>0$ and a set
$B\subseteq\mathbb{R}^{2}_{+}$ (having positive Lebesgue measure) such that
the function

$f(x)=\begin{cases}ax,&0\leq x\leq b\\\ -dx+(a+d)b,&b\leq x\leq b+c\\\
rx-(r+d)(b+c)+(a+d)b,&x\geq b+c\end{cases}$

satisfies

${\rm I\kern-1.60004ptE}\bigl{(}f(x+X)-f(X)\bigr{)}\geq f(x),\quad\forall
x\in\mathbb{R}_{+},$

but $f(x+y)\leq f(x)+f(y)$ for $(x,y)\in B$.

###### Proof.

Observe that there exists $N\in\mathbb{N}$ s.t. $b/N<c$. This implies that any
$x\in(0,\tfrac{b}{2N})$ and $y\in(b,b+\tfrac{b}{2N})$ satisfies $f(x+y)\leq
f(x)+f(y)$ (take $B=(0,\tfrac{b}{2N})\times(b,b+\tfrac{b}{2N})$). Now, let
$\boldsymbol{\theta}:=(a,b,c,d,r)$ and $H(x,\boldsymbol{\theta}):={\rm
I\kern-1.60004ptE}\bigl{(}f(x+X)-f(X)\bigr{)}-f(x)$. Since $f(0)=0$, we
necessarily have that $H(0,\boldsymbol{\theta})=0$; and thus showing that
$H^{\prime}(x,\boldsymbol{\theta})\geq 0$ holds for some $a,b,c,d,r>0$ is
sufficient. Therefore, we have

$H^{\prime}(x,\boldsymbol{\theta})=\int^{\infty}_{0}f^{\prime}(x+y)e^{-y}\textnormal{d}y-f^{\prime}(x),$

where

$f^{\prime}(x)=\begin{cases}a,&0<x<b\\\ -d,&b<x<b+c\\\ r&x>b+c.\end{cases}$

It is now required to find the range of the parameters for which
$H^{\prime}(x,\boldsymbol{\theta})\geq 0$ is satisfied. First, for
$x\in[0,b]$,

$\displaystyle H^{\prime}(x,\boldsymbol{\theta})$
$\displaystyle=\int^{b-x}_{0}ae^{-y}\textnormal{d}y-\int_{b-x}^{b-x+c}de^{-y}\textnormal{d}y+\int_{b-x+c}^{\infty}re^{-y}\textnormal{d}y-a$
$\displaystyle=(r+d)e^{-(b-x+c)}-(a+d)e^{-(b-x)},$

and hence $H^{\prime}(x,\boldsymbol{\theta})\geq 0$ yields that

$c\leq\log(\tfrac{r+d}{a+d}).$ (3.1)

Performing similar calculations for $x\in[b,b+c]$ and $x\in[b+c,\infty)$
yields no further information about the parameters. Then, since there always
exists some $N\in\mathbb{N}$ such that $b/N<c$ and $a,d,r>0$ s.t.
$c\leq\log(\tfrac{r+d}{a+d})$, these parameters satisfy the conditions of the
lemma. ∎