# Precise Object Placement Using Force-Torque Feedback

Osher Lerner, Zachary Tam, Michael Equi  Department of Electrical Engineering
and Computer Science
University of California at Berkeley
Berkeley, CA 94720, USA
{oshlern, ztam<EMAIL_ADDRESS>
###### Abstract

Precise object manipulation and placement is a common problem for household
robots, surgery robots, and robots working on in-situ construction. Prior work
using computer vision, depth sensors, and reinforcement learning lacks the
ability to reactively recover from planning errors, execution errors, or
sensor noise. This work introduces a method that uses force-torque sensing to
robustly place objects in stable poses, even in adversarial environments. On
46 trials, our method finds success rates of 100% for basic stacking, and 17%
for cases requiring adjustment.

##  I INTRODUCTION

Stable object stacking is a common task for home robots, industrial/warehouse
robots, and exploring robots. For example, robots may be required to stack
dishes in a sink, cabinet, or dishwasher, or stably stack tools and parts on a
shelf. Furthermore, when exploring cities after a natural disaster, or when
exploring other planets, the ability to construct impromptu shelter from
provided rubble is critical.

When humans stack objects, the last mile of placement is guided primarily by
fine-grain force sensing and reactive control. This is especially apparent
when sight is removed or when trying to precariously balance complex objects.
In such actions, humans empirically tend to lightly release the grasped object
in order to sense the small net object wrench and compensate accordingly.

In this work, we use an ATI Axia80 wrist-mounted force-torque sensor to
emulate this behavior. We are able to sense when the robot’s net wrench from
the grasped object is zero, which by Newton’s Third Law means that the robot
is applying zero net wrench to the grasped object. At this point, the object
can safely be released. This known stable placement can be useful for rock
stacking, robust placement and stacking of delicate objects, and human-robot
handover. The work considers extensions of the grasp matrix/friction cone
paradigm of grasping presented in lectures in the latter half of the course.

Previous work in robotic stone stacking [4, 5, 6] uses vision based sensing to
find stable stacking placements, however it generally lacks the ability to
recover from visual input noise and robot kinematic error. This work seeks to
apply the force-torque feedback methods commonly used for peg insertion [1, 2]
to stone stacking in order to enable contact-based readjustment for stable
stacking.

In applying force-torque sensing to stone stacking, we attempt to minimize
accidental motion of the grasped object and the existing structure, while
simultaneously exploring the stacking pose space to find one with the optimal
wrench. We thus present a policy that estimates a gradient for pose
optimization to reduce external torques to zero, i.e., to find a stable
stacking pose.

##  II RELATED WORK

###  II-A Force-Torque Feedback

Wrist-mounted force torque sensors such as the ATI Axia80 are very common for
high precision and contact-critical tasks such as peg insertion and deburring
[1], [2], [3]. [1] uses impedance control and system force/torque modeling to
determine and zero out peg reaction forces. [2] extends prior work to a
deformable pegboard; because of the nonlinear dynamics, the authors use deep
reinforcement learning to model the problem. [3] uses fuzzy PI control in
order to maintain strong contact with the desired part for deburring. For peg
insertion tasks, the authors of [1] and [2] consider only the last mile of
contact adjustments, starting from a vision-based initial guess. In this work,
we take a similar approach of starting from a good initial contact and making
force-torque based adjustments in order to reach our goal. However, we apply
these concepts to the problem of stacking, where the goal is less well-defined
and it is important to minimize disturbances to the environment.

###  II-B Stone Stacking

Prior work on stone stacking uses various types of vision-based sensing to
find optimal stacks. [4] uses an objects with known 3D models and then
simulates possible stacks to find a stable one to execute. [5] uses a
Q-learning approach for stone stacking where object models are known and the
problem is reduced to 2D by using simple extruded polygons as stones. It tests
with physical run-throughs of simulated rollouts. [6] uses deep reinforcement
learning to stack 3D objects of varying irregularity in simulation from a
ground-truth depth map of the bottom of the held object and the top of the
existing structure. [7] uses vision sensing and sim-to-real transfer of a
reinforcement learning agent to form especially challenging stacks of known
polyhedra and cylinders. This paper extends these works by applying force-
torque sensing to accommodate for noise and inaccuracies in vision sensing and
object modeling.

###  II-C Alternative Stacking

We also note prior use of a novel end effector for grasping and stacking
objects. [8] includes a demonstration of a novel “bubble” gripper which senses
finger normal and shear forces in order to detect grasped object geometry and
external forces. This work demonstrates similar flat object stacking using a
simpler and more robust parallel-jaw gripper.

##  III METHODS

Figure 1: Left: The force diagram of a grasp before placement (normal forces
at the jaws omitted for clarity). Right: The new forces on the rock during a
place attempt, and the resulting torque.

###  III-A Pipeline

We devise an algorithm to pick, place, and iteratively adjust placement until
the rock can be released stably. We first implement a standard pick-and-place
algorithm to pick a rock from the ground and form an initial guess for
placement. The end effector is quasi-statically lowered onto this pose until
we detect significant resistance from the force-torque sensor. After halting
in this configuration, we use the force-torque readings to calculate the
approximate point of contact and its surface normal. We then slightly raise
the rock to relieve the forces and adjust the pose to place the center of mass
of the rock above the previous point of contact, and optionally further
navigate along the surface towards a flatter region. This process repeats
until the torques detected at the time of contact are sufficiently small, and
we can release the rock in a stable configuration.

This algorithm is orchestrated through ROS. The robotic arm receives position
and velocity commands in the world frame, and converts them into a joint-space
control loop, respecting physical and safety constraints on position, speed,
and acceleration. We use a node for computer vision to estimate the position
of rocks on the ground and an initial guess for placement on the tower, but
currently spoof the outputs with an approximate oracle. These outputs are
varied adversarially to measure the performance of the rest of the pipeline.
The force torque sensor node constantly publishes readings, and internally
handles calibration and reference frame conversions. Finally the gripper node
takes calibration, open, and close commands.

###  III-B Physics Model

For the workspace frame, we choose the robot’s base frame, labeled $B$. We
denote the wrist frame by $W$, the gripper-tip frame by $T$, the rock’s COM
(center of mass) frame by $R$, and the contact point frame as $C$.

The rock’s true COM is unknown, so we approximate it to be directly at the
center of the parallel jaws: $g_{RT}=I$. Given a better estimate of the COM,
we’d simply have to adjust our measurements by an additional frame conversion.
(One possible extension of this project would be to use force-torque readings
to estimate the COM of the object.)

As we lower the rock, it makes contact with the tower. We label the point of
contact (relative to the rock’s COM) by $\vec{r}$. At the point of contact,
the tower exerts a normal force $\vec{F}_{N}$ on the rock in the direction of
the surface normal (Fig. 2). This results in a torque on the rock’s COM given
by $\vec{\tau}=\vec{r}\times\vec{F}_{N}$.

Figure 2: Left: The simplified force diagram of an attempted rock placement,
used to compute the normal force at rest. Middle: The torque induced on the
rock COM, which we reverse engineer to find the contact point r from a
measured torque. Right: The normal force is along the surface normal of the
tower. Since the tower is likely convex around the point of contact, this can
be used to find the direction towards a flat peak.

In general, a wrench in frame $A$ can be specified as a vector concatenation
of torque and force $\mathcal{F}_{A}=\begin{bmatrix}\vec{\tau}\\\
\vec{F}\end{bmatrix}$. Assuming a rigid body connects $A$ to a frame $B$
(related by $g_{AB}\in SE(3)$), the wrench felt in frame $B$ can be computed
as [9]

$\mathcal{F}_{B}=\text{Ad}_{g_{AB}}^{T}\mathcal{F}_{A}$ (1)

As long as our rigidity assumptions hold (e.g. no slipping), the wrench on the
rock is sensed as a force and torque on the wrist by our sensor. We convert
this reading to the rock COM frame using (Eq. 1) and a pre-measured
transformation between the sensor and the gripper tip frames.

We wish to reconstruct the position of the contact point and orientation of
the surface from our sensor readings. Even though we derived how to compute
the wrench felt by the rock, this wrench is a product of 2 unknowns: the
contact point $\vec{r}$ and the normal force $\vec{F}_{N}$. To solve this, we
leverage our quasi-static assumption to claim the net force on the rock is
$\approx 0$ when we measure.

The forces acting on the rock are gravity, the push of the gripper, and the
normal force from the tower (see Fig. 2). The gravitational force we can
measure while holding the rock in a hovering position. The push of the gripper
we computed by converting the wrench from the force torque sensor. This means
we can derive the normal force from

$\vec{F}_{net}=\vec{F}_{g}+\vec{F}_{push}+\vec{F}_{N}=0.$ (2)

By computing $\vec{F}_{N}$, we know the direction of the surface normal of the
tower at the contact point. Now we seek $\vec{r}$, knowing
$\vec{\tau}=\vec{r}\times\vec{F}_{N}$. We can recover

$\vec{r}_{T}=\frac{1}{|\vec{F}_{N}|^{2}}\vec{F}_{N}\times\tau$ (3)

where $\vec{r}_{T}$ is the component of $\vec{r}$ within the tangent plane of
the tower surface. The component normal to the surface cannot be deduced it
from measuring forces and torques. Fortunately, for the purpose of selecting a
new placement we need only shift our pose along the tower surface tangent
plane, since any other pose would either not make contact with the tower or be
inside it.

###  III-C Surface Normal Optimization

The above approach for computing the displacement to the point of contact is
used in our optimization algorithm to select the next placement such that the
rock center of mass is directly over the previous contact point. For complex
objects, finding a balanced pose requires not only centering the COM above the
contact point, but locating a good contact point at which the contact surface
is flat.

To this end, we utilize our calculation of the surface normal. However,
finding a flat contact point on the tower fundamentally involves higher order
information about the surface: curvature describes the rate of change of the
surface normal. Since we gain no curvature information from our measurement
process, we must look to some useful assumptions. The point of contact is
almost certainly extremity, since it is the first point on the surface of the
tower we reach on a linear approach, particularly since we take a vertical
approach and the objects we handle are wide rocks. Around extremities, the
surface is convex, meaning it curves away in every direction. Since we seek
the direction along the surface where the normal vector is closest to
vertically upwards, the best guess is therefore the direction within the
tangent plane that has the largest z-component. Mathematically, using the
normal vector $\hat{n}=\vec{F}_{N}/|\vec{F}_{N}|$ this direction is given as

$\hat{d}=-\hat{n}\times(\hat{n}\times\hat{z})$ (4)

We can then simply add a small step in the direction of $\hat{d}$ to
$\vec{r}_{T}$ to determine our proposed shift towards a more stable
configuration.

###  III-D Calibration

Calibration is thoroughly critical to the success of this approach. We go
about calibration in two phases. The first measures the forces and torques on
the wrist when the object is grasped and suspended statically by the robot.
The second measures the forces and torques on the wrist when putting the
object onto a flat surface. Taking a baseline measurement of the forces and
torques when pushing the object onto a flat surface, such as the location
where the object was picked up, gives us a reference point for when the object
is in a stable position. We can then subtract this reference from the forces
and torques when stacking to determine when the object is close enough to the
stable reference to be released without falling. The primary purpose of the
first measurement is to find the weight of the object / gripper combination
which is used to separate out the normal force as a result of gravity and the
force with which the robot is pushing down on the object.

When taking each force-torque measurement we wait for 0.1 seconds after the
motion of the arm has stopped and take an additional 0.5 seconds of samples
which are then averaged to devise our final force-torque reading at that point
in time. Without this step, our readings would be overwhelmed by noise and the
25Hz read-rate of the sensor. This is to help improve the consistency of the
measurements and thus result in faster convergence to stable placement
positions. (See Fig. 3)

Figure 3: Norms of force and torque readings upon contact during vertical
descent. Left: Centered placement. Right: Off-center placement.

Rigidity is central to our mathematical derivations, yet slippage is
commonplace between the object and the gripper as well as the object and the
tower. We designed our algorithm for consistency since our physical
assumptions were often violated. By approaching placement vertically each
time, we avoid slippage along the surface of contact. This also serves to zero
out deviations from our model due to misalignment of normal forces, gravity,
and the arm’s pushing forces. Due to imprecise grasping of the objects and the
nonlinear fashion in which our gripper pose depends on the object it is
gripping, our dynamical models failed to accurately compute terms and convert
between frames, but these complications are zeroed out when everything aligns
with the direction of gravity.

##  IV RESULTS

###  IV-A Physical Setup

In this work, we use a 2014 UR5 CB3 Cobot equipped with an ATI Axia80 force-
torque sensor and a Robotiq 2F-140 parallel-jaw gripper. Since this paper is
concerned primarily with the placement section of pick and place, we define a
preset picking process. As a stand-in for a vision-based object detection and
grasping algorithm, a human operator places objects in a specified pickup
location on flat ground. The gripper is able to repeatedly collect objects
from this loading zone and place them in a specified initial stacking
position. The initial stacking position is a hard-coded stand-in for an
optimization warm-start we could otherwise source from an algorithm like the
one presented in [6]. The existing stack structure is aligned with this
initial position before experiments begin.

Due to the difficulties we faced in the later parts of this work, we stack
with largely flat objects such as shuffleboard pucks (which have flat bottoms
but only undulating semi-flat tops) and flat squares of 1/4 inch fiberboard.

In these experiments, our inputs were the initial position of the existing
stack and the selected object to place, while our output was the success of an
attempted stack placement (a binary variable). We hypothesized that more
complex objects and existing stacks would result in lower likelihood of a
successful stack and more update steps before a stable pose was found.

###  IV-B Stacking Experiments

Using the setup described above, we performed 46 documented stacking tests. On
the 16 trials where no update steps are required (i.e., the initial stacking
pose corresponds to a stable placement), our policy detects this with 100%
accuracy and precision. This is successful even though the algorithm has no
information on how high the existing stack is (other than an upper bound). We
were able to stack at least 6 objects before we ran out of objects to stack in
the lab space (Fig. 4).

Figure 4: A six-object stack constructed from the ground-up by our policy.

However, we find only a 17% success rate on more difficult experiments.

On trials where the placement surface is non-horizontal, the algorithm will
not succeed unless we allow the robot to adjust the orientation of the end-
effector (not accounted for in this work). Instead, it makes the inaccurate
implicit assumption that the ramp leads to a stable placement at the top, thus
never finding a stable stacking pose on the ramp.

Even with near-horizontal stable stacking poses, when multiple steps are
required for stacking we also find limited success (18%). In practice, we find
the gradient update steps to be inconsistent at moving the object toward a
stable pose, and we find that the metric for finding a stable placement pose
is often misaligned with the real world.

###  IV-C Algorithm Component Experiments and Analysis

Due to the challenges we faced in actual stacking experiments, we performed
other experiments and analysis with smaller components of our algorithm. These
results serve to demonstrate the individual success and theoretical
feasibility of these components, and help to narrow down the source of any
remaining error.

In early development before the ROS architecture properly communicated, we
conducted experiments to manually analyze and visualize the force-torque
readings from place actions. This procedure proved insightful later and
throughout the project to diagnose the algorithm’s performance. We list
example data from three attempted placements of a held puck unto a puck on the
ground. The readings are calibrated with the object suspended in the air, and
processed to eliminate remaining offsets.

Figure 5: Left: Experimental readings of characteristic states from 3 place
maneuvers. Upon first contact, noise still overwhelms torque readings. After
pressing with 10N, the direction of the stable site can be inferred from the
data. Right: Plotted center of mass for each candidate pose and the
algorithm’s suggested shift. ”Centered” in green, ”Y-shift” in red, ”X-shift”
in blue, and true tower center in black.

In the first experiment, we spoof a ”guessed” placement pose from vision to be
quite close to the true center of the puck below it. The second is offset
mostly in the y direction, and the third is offset mostly in the x direction.
The data from characteristic states during placement is listed in (Fig 5), and
the initial placements relative to the tower center (0,0) are plotted. Once
sufficiently pressed against the tower, we plug the force torque readings into
our algorithm to propose candidate shifts. The proposed shifts are remarkably
successful, pointing to locations within 1cm of the optimal placement in just
one iteration.

However, we notice the proposed shift for the already nearly centered
placement is an arbitrary direction, leading to a worse guess for the next
attempt. We observe this failure mode consistently throughout our experiments:
update steps are accurately computed for egregiously offset placements, but as
we near a stable configuration the algorithm fails to converge and proposes
further poses in random directions.

When the pipeline failed to place objects, we would create manual contact
forces to verify the accuracy of our readings, calibration, and update steps.
This process is done by simply applying a strong pressure to the object with
the experimenter’s finger while the robotic arm lowers it towards the tower.
We were able to get clear readings with up to 30Nm of torque, and the proposed
shift directions passed this sanity check 100% of the time.

###  IV-D Implications

Given the independent success of our placement algorithm and our update steps,
both in simulation and in sanity checking, we are able to narrow down the
source of the error when we combine the subcomponents of our method.

While the finger pressing sanity check demonstrates the success of our
algorithm for high forces, it says little about how well it scales down (i.e.,
for the very small torques when the current pose is only marginally off from a
stable pose). Furthermore, stacking tests where update steps were not required
seem to show that the algorithm is able to detect stable poses. Therefore, we
suspect a sensitivity limit to our algorithm, i.e., it cannot detect the
direction of small torques with enough accuracy to get the object into a
stable pose while simultaneously being able to distinguish a stable pose from
one that is unstable but nearby. This explains the observed pattern of
behavior of accurate adjustments from extremely unstable placements, then
incorrect adjustments from partially stable poses.

##  V CONCLUSION AND FUTURE WORK

In this work, we were able to use solely force-torque sensing to consistently
stack reasonably flat objects in a given location and to occasionally find
stable stacks from an unstable initial pose.

When working on this project, we encountered two major difficulties. Firstly,
we were trying to implement novel software on three new independent pieces of
hardware. Ultimately, setup took a considerable amount of time, as did tuning
the force-torque sensor (Section III-D). Once we had set up the hardware, we
encountered significant challenges with getting proper force-torque readings
in order to find accurate gradient steps (see Section IV).

The first steps for future work are to incorporate better methods for
detecting small offsets from stable poses; this is where we found difficulty
in this work. In order to leverage our successful approach into consistent
convergence, we must modify our calculations for nearly-balanced poses.

One possibility is to slowly release the objects in order to predict how they
might behave if fully released. We can thus determine whether to release or
where to move the object next.

Once we are able to consistently stack flat rocks, we plan to stack real rocks
which are bulky, irregular, and rough. This comes with nonconvexity, multiple-
contact points, and ill-posed optimization. For these cases we can utilize our
surface normal optimization edition to our procedure. Although our derivation
assumed many good conditions, we believe it will succeed on most common
objects, particularly due to our resistance threshold and the nature of
placing.

For objects that slip into unstable configurations, or generally rotated pick
and place actions, we need to account for more than just vertical motion. This
requires consideration of the gripper model center of mass and the wrench it
applies to the force-torque sensor. It also requires optimizing over a higher
dimensional space, where we must consider out-of-plane torques and changes in
gripper orientation as optimization steps.

For multi-rock sequential stacking, we can record the final position of the
center of mass of each placed rock and use it as our initial guess for the
next rock’s placement. When our control loop fails to converge, we place the
rock on the ground and grasp the last rock on the tower using its stored pose.
We remove it as well, and try a random new sequence of rock stacks, until
those are exhausted and we backtrack further down the tower. This search
doubles as a reordering of rocks as well as just random repositioning due to
the noise in our stacking algorithm.

A major challenge we encountered was with the stickiness of the gripper
fingers. In order to press against the environment without losing grip of the
target object, high friction is required for the gripper fingers. However,
this comes with stickiness, which tends to cause the grasped object to move
upon release (i.e., it gets stuck to the gripper fingers). This problem is
especially apparent with the slick shuffleboard pucks.

A promising future approach to tackling the complex nonlinear nature of this
optimization procedure is to leverage deep learning. Some assumptions have to
be made to make progress in the pick-and-place task, and deep learning enables
us to implicitly learn the relevant assumptions from patterns in data. This
can be beneficial for both general agents with hard-to-specify assumptions and
specialized datasets for particular objects and environments. Today’s cutting
edge research [5] uses reinforcement learning, which our method could utilize
to form strategies for sequentially sampled placements; incorporating an RL
agent could enable interaction with the environment with the intent of
information gathering to quicken long-term success, as opposed to only
attempting successful placement with the present action. A deep model could
also potentially hybridize the tasks of vision and tactile sensing, leveraging
the relationship between the two more than a separated pipeline could.

Alternatively, a 0th order sampling of higher order information is possible:
by sampling readings at a local array of points, we can interpolate gradients
towards flat and balanced positions.

##  VI ACKNOWLEDGMENTS

The authors gratefully acknowledge the contribution of National Research
Organization and reviewers’ comments. We would also like to thank Professor
Ruzena Bajcsy, Isabella Huang, and Ellis Ratner for generously allowing us to
borrow their equipment and laboratory space.

## References

  * [1] J. Luo, E. Solowjow, C. Wen, J.A. Ojea, A.M. Agogino, “Deep Reinforcement Learning for Robotic Assembly of Mixed Deformable and Rigid Objects”, in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Madrid, Spain, October 1-5, 2018.
  * [2] P. Zou, Q. Zhu, J. Wu, J. Jin, “An Approach for Peg-in-Hole Assembling Based on Force Feedback Control”, in IEEE Chinese Automation Conference, Hangzhou, China, November 22-24, 2019.
  * [3] J.N. Pires, J. Ramming, S. Rauch, R. Araujo, “Force/torque sensing applied to industrial robotic deburring”, in Sensor Review, Bradford, United Kingdom, 2002.
  * [4] F. Furrer, M. Wermelinger, H. Yoshida, F. Gramazio, M. Kohler, R. Siegwart, M. Hutter, “Autonomous robotic stone stacking with online next best object target pose planning”, in IEEE International Conference on Robotics and Automation (ICRA), Singapore, May 29-June 3, 2017.
  * [5] Y. Liu, S.M. Shamsi, L. Fang, C. Chen, N. Napp, “Deep Q-Learning for Dry Stacking Irregular Objects”, in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Madrid, Spain, October 1-5, 2018.
  * [6] A. Menezes, P. Vicente, A. Bernardino, R. Ventura, “From Rocks to Walls: a Model-free Reinforcement Learning Approach to Dry Stacking with Irregular Rocks”, in IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), Nashville, TN, June 19-25, 2021.
  * [7] A.X. Lee, C. Devin, Y. Zhou, T. Lampe, K. Bousmalis, J.T. Springenberg, A. Byravan, A. Abdolmaleki, N. Gileadi, D. Khosid, C. Fantacci, J.E. Chen, A. Raju, R. Jeong, M. Neunert, A. Laurens, S. Saliceti, F. Casarini, M. Riedmiller, R. Hadsell, F. Nori, in Conference on Robot Learning (CoRL), London, United Kingdom, November 8-11, 2021.
  * [8] N. Kuppuswamy, A. Alspach, A. Uttamchandani, S. Creasey, T. Ikeda, R. Tedrake, “Soft-bubble grippers for robust and perceptive manipulation”, in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Las Vegas, NV, October 25-29, 2021.
  * [9] Lynch, Kevin M, and Frank C Park. “Wrenches.” Modern Robotics, Northwestern University, 26 Aug. 2017, https://modernrobotics.northwestern.edu/nu-gm-book-resource/3-4-wrenches/.


monthdayyeardate[] ,

# On the Density arising from the Domain of Attraction between Sum and
Supremum: the $\alpha$-Sun operator

N.S. Witte  School of Mathematics and Statistics, Victoria University of
Wellington, PO Box 600 Wellington 6140, New Zealand<EMAIL_ADDRESS>and  P.E. Greenwood  Department of Statistics, University of British Columbia,
Vancouver BC V6T 1Z4, Canada<EMAIL_ADDRESS>
###### Abstract.

We explore the analytic properties of the density function
$h(x;\gamma,\alpha)$, $x\in(0,\infty)$, $\gamma>0$, $0<\alpha<1$ which arises
from the domain of attraction problem for a statistic interpolating between
the supremum and sum of random variables. The parameter $\alpha$ controls the
interpolation between these two cases, while $\gamma$ parametrises the type of
extreme value distribution from which the underlying random variables are
drawn from. For $\alpha=0$ the Fréchet density applies, whereas for $\alpha=1$
we identify a particular Fox H-function, which are a natural extension of
hypergeometric functions into the realm of fractional calculus. In contrast
for intermediate $\alpha$ an entirely new function appears, which is not one
of the extensions to the hypergeometric function considered to date. We derive
series, integral and continued fraction representations of this latter
function.

###### 2010 Mathematics Subject Classification:

60E07,60G70,30A80,33.20,33A35

##  1\. Our Motivation

Consider a sequence of independent random variables $X_{1},X_{2},\ldots$ with
a distribution $F\in D(\Phi_{\gamma})$, where $D(.)$ is the domain of
attraction of the distribution $\Phi_{\gamma}$ and this latter distribution is
the extreme value distribution
$\Phi_{\gamma}(x)=\exp(-x^{-\gamma})\mathbb{I}_{[0,\infty)}(x)$. Now define
the following sequence $Y_{0},Y_{1},Y_{2},\ldots$

(1.1) $Y_{n}=\max\\{Y_{n-1},\alpha Y_{n-1}+X_{n}\\},\quad n\geq 1,\quad
Y_{0}=X_{0}\in\mathbb{R}.$

For $\alpha=0$ we have the supremum,
$Y_{n}=\max\\{X_{0},X_{1},X_{2},\ldots,X_{n}\\}$, while for $\alpha=1$ we have
the sum, $Y_{n}=X_{0}+\sum^{n}_{j=1}X_{j}\mathbb{I}_{[0,\infty)}(X_{j})$. In
this way one can interpolate between the sum and the supremum of a sequence of
i.i.d. random variables via the parameter $\alpha$.

There is a compelling motivation for our model (1.1) and the associated
operator interpolating between sum and supremum. Many decisions which we take
are of the following form: Accept the status quo, in which case
$Y_{n}=Y_{n-1}$, or risk losing some portion of what we have, in the process
of taking some action, so that the value of our current state is reduced to a
portion $\alpha Y_{n-1}$, in which case we can increase our value by a
stochastic $X_{n}$ to obtain $Y_{n}=\alpha Y_{n-1}+X_{n}$. Suppose we are so
lucky or wise that we can evaluate $X_{n}$ before taking the decision whether
to act or not. Then we can maximize the outcome by formulating our decision as
per (1.1).

In their study of this model - the $\alpha$-sun process - Greenwood and
Hooghiemstra [18, Theorem 2, Eq. (2.4)] derived a linear, homogeneous integral
equation for the density

(1.2) $h(x)=\frac{\gamma}{x}\int_{0}^{x}du\;\frac{h(u)}{(x-\alpha
u)^{\gamma}},$

subject to the following conditions: $\alpha\in(0,1)$, $\gamma\in(0,\infty)$
and $h(x)$ is a real, normalised density on $x\in(0,\infty)$. No further
conditions can apply. We investigate the solutions $h(x)=h(x;\gamma,\alpha)$
to this integral equation here.

In a subsequent work Hooghiemstra and Greenwood [25] took up the related
problem of finding the domain of attraction for the other extreme value
distributions, namely $F\in D(\Psi_{\gamma})\cup D(\Lambda)$ where
$\Psi_{\gamma}(x)=\exp(-(-x)^{\gamma})\mathbb{I}_{(-\infty,0)}(x)+\mathbb{I}_{(0,\infty)}(x)$
and $\Lambda(x)=\exp(-e^{-x})$. For the latter case a solution for the density
was given there by Eq. (10) in Cor. 2 of §2,

(1.3)
$h(x;\alpha)=\frac{1-\alpha}{\Gamma((1-\alpha)^{-1})}\exp\left(-x-e^{-(1-\alpha)x}\right).$

In the former case one seeks solutions to an integral equation similar to
(1.2), as given by Eq. (9) in Cor. 1,

(1.4) $h(x)=\frac{\gamma}{|x|}\int_{x/\alpha}^{x}du\;|x-\alpha
u|^{\gamma}h(u),\quad x\in(-\infty,0).$

However no explicit solutions to either (1.2) nor (1.4) have been reported in
any studies following this line of enquiry.

Using different methods Schlather [48] has treated a related problem. The
Schlather paper says, at the very beginning, that the 1991 and 1997 Greenwood
and Hooghiemstra papers combine the CLT with extreme value theory. However
this is not the case. We take up this distinction in the discussion §5. See
the more recent thesis by Anja Janßen [27] continuing this other line of
enquiry.

The study of density solutions to a variant of (1.2), specialised to
$\alpha=1$ but possessing an additional parameter $m\in\mathbb{N}$, was
initiated in [49], then pursued in [39] and most recently by [28]. This
distinct integral equation takes the form

(1.5) $h(x)=\frac{\gamma}{x^{m}}\int_{0}^{x}du\;\frac{h(u)}{(x-u)^{\gamma}},$

and the corresponding random variables are called ”generalized stable” and our
case includes the classical stable case $m=1$. Such random variables are of
interest in the case $m=2$, $\gamma\in(0,1)$ which is especially investigated
in Section 3 of [49] and Section 7 of [39], because of its relevance to
particle transport on a one-dimensional lattice. The first work [49]
identifies the solutions in terms of Fox $H$-functions and shows that for all
$m\in\mathbb{N}$, $\gamma\in(1-m,1)$ there exists a density solution to (1.5)
having a convergent power series representation at infinity and a Fréchet-like
behaviour at zero. The second work [39] considers invariance under length-
biasing and shows that for all $m\in\mathbb{N}$, $\gamma\in(1-m,1)$ there
exists a unique density solution to (1.5), whose corresponding random variable
can be represented in the case $\gamma\in(1-m,2-m)$ as a finite independent
product involving the Gamma and positive stable random variables. We have
several common results with [28], which are given later, and this latter work
shows that there exists a density solution, which is then unique and can be
expressed in terms of the Beta distribution, if and only if $m>1-\gamma$.
These density solutions extend the class of generalised one-sided stable
distributions introduced in [49] and subsequently investigated in [39]. This
work studies various analytical aspects of these densities, and solves some
open problems about infinite divisibility formulated in [49].

From the analytical viewpoint there is another thread that our study connects
with, and this concerns the solutions of integro-differential equations or
more specifically differential-delay Volterra integral equations. For this to
become apparent let us consider the cumulative distribution functions

(1.6) $\displaystyle f(x):=\int_{0}^{x}du\;h(u),\quad
f(0)=0,f(\infty)=1,\quad\text{from the solution of \eqref{integral_eqn}},$
(1.7) $\displaystyle f(x):=\int_{x}^{\infty}h(-u),\quad
f(0)=1,f(\infty)=0,\quad\text{from the solution of \eqref{integral_eqn-
variant}}.$

An important point to note is that henceforth we will frame our discussion of
(1.4) for $x\in\mathbb{R}_{+}$ by defining another density $h_{<}(x):=h(-x)$.
Where there is no chance of confusion we will often drop the subscript $<$.
Thus our defining integral equation will become

(1.8) $h_{<}(x)=\frac{\gamma}{x}\int_{x}^{x/\alpha}du\;(x-\alpha
u)^{\gamma}h_{<}(u),\quad x\in(0,\infty).$

A simple exercise then recasts (1.2) into the retarded version of such
equations

(1.9) $f^{\prime}(x)-(1-\alpha)^{-\gamma}\gamma
x^{-\gamma-1}f(x)=-\alpha\gamma^{2}x^{-\gamma-1}\int_{0}^{1}dq\;(1-\alpha
q)^{-\gamma-1}f(qx),$

whereas (1.4) becomes the advanced version

(1.10) $f^{\prime}(x)+(1-\alpha)^{\gamma}\gamma
x^{\gamma-1}f(x)=\alpha\gamma^{2}x^{\gamma-1}\int_{1}^{\alpha^{-1}}dq\;(1-\alpha
q)^{\gamma-1}f(qx),$

where these are to solved as boundary value problems, as defined above. This
pair of equations is a specialised form of a general class of integro-
differential equations formulated by Iserles and Liu [26], although they did
not derive solutions of the general system nor our particular case. This class
also covers the Pantograph type equations [11], [3] and cell division and
growth models [21], [22]. An early study presaging these developments is one
on a stochastic absorption problem [16].

We will only treat the solutions of (1.2) or (1.9) in this work, and defer the
solutions of (1.8) or (1.10) to a subsequent study. Following this
introduction we find explicit solutions to the two special cases $\alpha=0$
and $\alpha=1$ in §2. The density functions we find here are either well-known
or fairly well studied extensions to the hypergeometric functions, and we
include these cases because they serve as ”bookends” to the novel aspects of
our results. After that we treat the general case $0<\alpha<1$ in §3, which
can’t be handled in the same way as the preceding section. Our main results
are given in Prop. 3.9. The density function arising here is somewhat more
exotic, one of a class of functions lying beyond the generalised
hypergeometric, Meijer-G or Fox H-functions. We give explicit series, integral
and continued fraction representations of the specific function arising in our
study. In §4 we evaluate and plot the densities for a range of
$\gamma\in\\{\tfrac{1}{4},\tfrac{1}{2},\tfrac{3}{4},1\\}$ and
$\alpha\in\\{\tfrac{1}{4},\tfrac{1}{2},\tfrac{3}{4},1\\}$. A number of
implications and conclusions are taken up in the Discussion §5.

##  2\. Special Cases $\alpha=0$ and $\alpha=1$

###  2.1. Case $\alpha=0$: the Fréchet and Weibull distributions

In this case (1.2) reduces to

(2.1) $h_{0}(x)=\frac{\gamma}{x^{\gamma+1}}\int_{0}^{x}du\;h_{0}(u).$

A simple differentiation of this yields the ordinary differential equation

$\frac{d}{dx}h_{0}+\left[\frac{\gamma+1}{x}-\frac{\gamma}{x^{\gamma+1}}\right]h_{0}=0,$

with the solution

$h_{0}=Cx^{-\gamma-1}\exp(-x^{-\gamma}),$

and the normalisation gives $C=\gamma$. Thus

(2.2) $h_{0}(x)=\gamma
x^{-\gamma-1}\exp\left[-x^{-\gamma}\right]=\frac{d}{dx}\exp\left[-x^{-\gamma}\right].$

See [15], [14], [20], [47].

It will be of interest to compute the Mellin transform of $h_{0}(x)$ as given
by (2.2), which we denote by $H_{0}(s)$ for $\Re(s)<1+\gamma$

(2.3) $H_{0}(s)=\Gamma\left(\frac{1+\gamma-s}{\gamma}\right).$

Thus we can construct $h_{0}(x)$ from the inverse Mellin transform

$h_{0}(x)=\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}ds\;x^{-s}\Gamma\left(\frac{1+\gamma-s}{\gamma}\right),\quad
c<1+\gamma.$

###  2.2. Case $\alpha=1$

In the retarded case (1.2) reduces to a linear, homogeneous and singular
integral equation

(2.4)
$h_{1}(x)=\frac{\gamma}{x}\int_{0}^{x}du\;\frac{h_{1}(u)}{(x-u)^{\gamma}}.$

######  Proposition 2.1.

Let $0<\Re(\gamma)<1$. Then the solution $h_{1}(x)$ has the integral
representation

(2.5)
$h_{1}(x)=\frac{1}{\pi}\int_{0}^{\infty}dr\;\exp\left[-xr-\cos(\pi\gamma)\Gamma(1-\gamma)r^{\gamma}\right]\sin\left(\frac{\pi}{\Gamma(\gamma)}r^{\gamma}\right).$

###### Proof.

Being a singular integral equation with difference kernel, (2.4) can be
treated with a Laplace transform. For $\Re(p)>0$ let

$L_{1}(p):=\int_{0}^{\infty}dx\;e^{-px}h_{1}(x).$

Then the Laplace transform of (2.4) gives us

$\displaystyle-\gamma^{-1}\frac{d}{dp}L_{1}(p)$
$\displaystyle=\int_{0}^{\infty}dx\;\int_{0}^{x}du\;e^{-px}h_{1}(u)(x-u)^{-\gamma},$
$\displaystyle=\int_{0}^{\infty}du\;h_{1}(u)\int_{u}^{\infty}dx\;e^{-px}(x-u)^{-\gamma},$
$\displaystyle=\Gamma(1-\gamma)p^{\gamma-1}L_{1}(p).$

This latter equation is easily solved for $L_{1}(p)$ and upon the inversion of
the transform yields

(2.6)
$h_{1}(x)=\frac{C}{2\pi}\int_{c-i\infty}^{c+i\infty}dp\;e^{px-\Gamma(1-\gamma)p^{\gamma}},\quad
c>0,$

where $C$ is a normalisation constant. The integrand of the inversion formula
has no non-analytic features other than a branch cut in the $p$-plane on
$[0,-\infty)$ and it is possible to deform the vertical contour to a Hankel
loop $(-\infty,0,-\infty)$ either side of the cut. Subsequently, in the second
proof of Prop. 3.4, we will see that $C=1$. This then gives the stated result,
(2.5). ∎

######  Remark 2.1.

When $\alpha=1$ we are inevitably lead to the additional constraint of
$\Re(\gamma)<1$. This is clear from the strength of the singularity at the
end-point of the integration range in (2.4) and Gamma function
$\Gamma(1-\gamma)$ in (2.6)

######  Proposition 2.2.

For $\Re(\gamma)<1$ the density $h_{1}(x)$ has a convergent large-$x$
expansion

(2.7) $h_{1}(x)=\frac{1}{\pi x}\sum_{n=1}^{\infty}(-1)^{n-1}\sin(\pi
n\gamma)\frac{\Gamma(1+n\gamma)}{n!}\left(\frac{\Gamma(1-\gamma)}{x^{\gamma}}\right)^{n}.$

###### Proof.

Rewriting the integrand of (2.5) in terms of exponentials instead of
trigonometric functions and then performing the Taylor expansion of
$\exp\left(-e^{\pm\pi i\gamma}\Gamma(1-\gamma)r^{\gamma}\right)$ about
$r^{\gamma}=0$ the resulting integrals can be done, after recognising the sum
and integral are uniformly and absolutely convergent. ∎

######  Remark 2.2.

Formula (2.7) can be simplified when $\gamma\in\mathbb{Q}$. If $\gamma$ is
rational then Gauss’s multiplication formula for the Gamma function [45, Eq.
(5.5.6)] can be employed to write the sum in a hypergeometric form. For
example one has

(2.8) $\displaystyle h_{1}(x;\tfrac{1}{2})$
$\displaystyle=\frac{1}{2x^{3/2}}e^{-\pi/4x},$ (2.9) $\displaystyle
h_{1}(x;\tfrac{1}{3})$
$\displaystyle=\frac{\Gamma\left(\frac{2}{3}\right)}{\sqrt[3]{3}x^{4/3}}\text{Ai}\left(3^{-1/3}\Gamma(\tfrac{2}{3})x^{-1/3}\right),$
(2.10) $\displaystyle h_{1}(x;\tfrac{2}{3})$ $\displaystyle=\frac{4\pi
e^{-2\Gamma\left(\frac{4}{3}\right)^{3}/x^{2}}}{3^{5/6}\Gamma\left(-\frac{1}{3}\right)x^{7/3}}$
$\displaystyle\times\left[3^{2/3}x^{2/3}\text{Ai}^{\prime}\left(3^{2/3}\Gamma\left(\tfrac{4}{3}\right)^{2}x^{-4/3}\right)-3\Gamma\left(\tfrac{4}{3}\right)\text{Ai}\left(3^{2/3}\Gamma\left(\tfrac{4}{3}\right)^{2}x^{-4/3}\right)\right],$

using the standard definitions of Airy function [45, Chapter 9].

The following result for general $\gamma$ appears to be simplest
identification amongst the extensions of hypergeometric function families.

######  Proposition 2.3 (§3.1.3, [28]).

The density $h_{1}(x;\gamma)$ is the Fox $H$-function

(2.11)
$h_{1}(x;\gamma)=\gamma^{-1}\left(\Gamma(1-\gamma)\right)^{-1/\gamma}{\rm
H}^{0,1}_{1,1}\left(\frac{x}{\left(\Gamma(1-\gamma)\right)^{1/\gamma}};\begin{array}[]{c}(1-\gamma^{-1},\gamma^{-1})\\\
(0,1)\end{array}\right),$

as per the definition of [30][Eq. (51), pg. 316, Chap. 7].

###### Proof.

This follows from the Mellin-Barnes integral representation of the Fox
$H$-function, see (3.11). This also agrees with the formulae in §3.1.3 of
[28], with $m=n=1$ and the scale factor
$c=\left(\gamma\Gamma(1-\gamma)\right)^{-1/\gamma}$. ∎

Figure 1. Densities $h_{1}(x;\gamma)$ versus $x$ for
$\gamma=\\{\frac{1}{4},\frac{1}{3},\frac{1}{2},\frac{2}{3},\frac{3}{4}\\}$.

##  3\. General $\alpha$ Case

The kernel in (1.2) is not separable nor of convolution type unless $\alpha=1$
so we need to employ more powerful methods in the generic $\alpha$ case, than
those utilised in §2. In the general case the Mellin transform of the density
$h(x)$ is defined by

(3.1) $H(s):=\int_{0}^{\infty}dx\;x^{s-1}h(x),$

for a suitable vertical strip in the $s$-plane, to be clarified below, but
containing the line $\Re(s)=1$. This constitutes our normalisation $H(1)=1$.

######  Proposition 3.1.

For $0<\alpha<1$, $\gamma>0$ and $\Re(s)<1+\gamma$ the retarded Mellin
transform $H(s;\gamma,\alpha)$ satisfies the linear, homogeneous functional
equation

(3.2)
$H(s)=\frac{\gamma}{1+\gamma-s}{}_{2}F_{1}(\gamma,1+\gamma-s;2+\gamma-s;\alpha)H(s-\gamma),$

where ${}_{2}F_{1}(a,b;c;z)$ is the Gauss hypergeometric function [45, §15].

###### Proof.

Taking the Mellin transform of both sides of (1.2) we compute

$\displaystyle\int_{0}^{\infty}dx\;x^{s-1}h(x)$
$\displaystyle=\gamma\int_{0}^{\infty}dx\;x^{s-2}\int_{0}^{x}du\;\frac{h(u)}{(x-\alpha
u)^{\gamma}},$
$\displaystyle=\gamma\int_{0}^{\infty}du\;h(u)\int_{u}^{\infty}dx\;\frac{x^{s-2}}{(x-\alpha
u)^{\gamma}},$
$\displaystyle=\gamma\int_{0}^{\infty}du\;h(u)u^{s-1-\gamma}\int_{1}^{\infty}dx\;\frac{x^{s-2}}{(x-\alpha)^{\gamma}}.$

The nested integral exists only for $\Re(s)<1+\gamma$ and the integrand
permits an expansion about $\alpha=0$, which can be integrated term-by-term to
yield

$\int_{1}^{\infty}dx\;\frac{x^{s-2}}{(x-\alpha)^{\gamma}}=\frac{1}{1+\gamma-s}\sum^{\infty}_{k=0}\frac{(\gamma)_{k}(1+\gamma-s)_{k}}{k!(2+\gamma-s)_{k}}\alpha^{k}.$

This allows us to make the identification with the stated Gauss hypergeometric
function and (3.2).

From its definition (3.1) $H(s)$ exists for $s$ in the strip $a<\Re(s)<b$, for
some $a<b$. Due to the existence of the norm we have $a<1<b$. In addition
$H(s-\gamma)$ exists for $a<\Re(s)-\gamma<b$, or $a+\gamma<\Re(s)<b+\gamma$.
There must be a common overlap of these two intervals so we conclude
$a+\gamma<b$. Exploiting all of the above inequalities we find that
$b<1+\gamma$ and $1-\gamma<a$ and so in summary we deduce
$1-\gamma<a<1<b<1+\gamma$. ∎

###  3.1. On the nature of the
${}_{2}F_{1}(\gamma,1+\gamma-s;2+\gamma-s;\alpha)$

Firstly with regard to the analytical character of this hypergeometric
function with respect to $s$ we can make the following immediate observations:

  1. (a)

For $0<|\alpha|<1$ the ${}_{2}F_{1}$ is convergent for all $\Re(s)<1+\gamma$.
The ${}_{2}F_{1}$ has simple poles on the right at $s=1+\gamma+k$ for
$k\in\mathbb{Z}_{\geq 0}$.

  2. (b)

For $|\alpha|=1$ we require in addition that $\Re(\gamma)<1$. The poles are
the same as for $|\alpha|<1$,

  3. (c)

Our ${}_{2}F_{1}$ function is one of a pair of solutions to the hypergeometric
differential equation about $\alpha=0$ that is bounded as $s\to-\infty$ (by
$(1-\alpha)^{-\gamma}$). This distinguishes itself from the other member of
the pair which has the simple algebraic form $\alpha^{s-\gamma-1}$, and which
diverges as $s\to-\infty$.

This hypergeometric function can linearly mapped into other forms via the
Kummer or linear fractional transformations of the $\alpha$-plane, and we
record three versions [45, 15.8.1] of these

$\displaystyle{}_{2}F_{1}(\gamma,1+\gamma-s;2+\gamma-s;\alpha)$
$\displaystyle=(1-\alpha)^{-\gamma}{}_{2}F_{1}(\gamma,1;2+\gamma-s;-\tfrac{\alpha}{1-\alpha}),$
$\displaystyle=(1-\alpha)^{s-\gamma-1}{}_{2}F_{1}(2-s,1+\gamma-s;2+\gamma-s;-\tfrac{\alpha}{1-\alpha}),$
$\displaystyle=(1-\alpha)^{1-\gamma}{}_{2}F_{1}(2-s,1;2+\gamma-s;\alpha).$

The hypergeometric function is a special case of the most general, because
$c=b+1$ and only contains one general argument and two general parameters.
Thus it can be identified in a number of ways. In the first instance it is an
incomplete Beta function, see [45, 8.17.7],

${}_{2}F_{1}(\gamma,1+\gamma-s;2+\gamma-s;\alpha)=(1+\gamma-s)\alpha^{s-\gamma-1}B_{\alpha}(1+\gamma-s,1-\gamma).$

In the second instance it is a generalisation of the Hurwitz-Lerch zeta
function, namely the higher-order Apostol-Lerch function defined by the
integral [33], [34], [24]

$\zeta(z,\sigma,v,\lambda)=\frac{1}{\Gamma(\sigma)}\int_{0}^{\infty}dx\;\frac{x^{\sigma-1}e^{-vx}}{\left(1-ze^{-x}\right)^{\lambda}},$

for $|z|\leq 1$, $\Re(v)>0$, $\Re(\sigma)>1$ and $\Re(\lambda)>0$. Thus

${}_{2}F_{1}(\gamma,1+\gamma-s;2+\gamma-s;\alpha)=(1+\gamma-s)\zeta(\alpha,1,1+\gamma-s,\gamma).$

A useful alternative integral representation of the ${}_{2}F_{1}$ is

(3.3)
${}_{2}F_{1}(\gamma,1;2+\gamma-s;\tfrac{\alpha}{\alpha-1})=1-\alpha\gamma(1-\alpha)^{\gamma}\int_{0}^{1}dt\;t^{1+\gamma-s}(1-\alpha
t)^{-\gamma-1}.$

###  3.2. Inequalities

We will require bounds on the density in our subsequent analysis, which can be
simply derived.

######  Proposition 3.2.

For $0<x<\infty$, $0<\alpha<1$ and $\gamma>0$ then retarded density
$h(x;\gamma,\alpha)$ is a real, positive density and satisfies the bounds

(3.4) $\gamma
x^{-\gamma-1}\exp\left[-(1-\alpha)^{-\gamma}x^{-\gamma}\right]\leq
h(x)\leq\gamma(1-\alpha)^{-\gamma}x^{-\gamma-1}\exp\left[-x^{-\gamma}\right].$

###### Proof.

We start by noting the simple inequality $x^{-\gamma}\leq|x-\alpha
u|^{-\gamma}\leq(1-\alpha)^{-\gamma}x^{-\gamma}$ for $0\leq u\leq x$.
Furthermore define the cumulative distribution function
$f(x):=\int_{0}^{x}du\;h(u)<1$ in the standard way. Employing this inequality
in (1.2) we have

$\gamma x^{-\gamma-1}f(x)\leq
h(x)=f^{\prime}(x)\leq\gamma(1-\alpha)^{-\gamma}x^{-\gamma-1}f(x).$

Integrating this from the reference point $x_{0}=\infty$ ($f(\infty)=1$) to
$x<\infty$ without crossing any standard placement of the branch cut, e.g.
$[0,-\infty)$, we deduce a reversed inequality

$\exp\left[-(1-\alpha)^{-\gamma}x^{-\gamma}\right]\leq
f(x)\leq\exp\left[-x^{-\gamma}\right].$

Combining these two inequalities we conclude (3.4). ∎

######  Corollary 3.1.

Let $\Re(s)=\sigma$. For $\sigma<1+\gamma$ the Mellin transform
$H(s;\gamma,\alpha)$ satisfies the bound

(3.5)
$|H(s)|\leq(1-\alpha)^{-\gamma}\Gamma\left(\frac{1+\gamma-\sigma}{\gamma}\right).$

###### Proof.

This follows immediately from (3.4). It should be noted that the stronger
result $|H(s)|\leq 1$ holds when $\Re(s)=1$. ∎

If one were to define a new transform $K(s)$ by

(3.6) $K(s):=\frac{(1-\alpha)^{s}}{\Gamma(1+\gamma^{-1}-\gamma^{-1}s)}H(s),$

then the functional equation becomes

(3.7)
$K(s)={}_{2}F_{1}(\gamma,1;2+\gamma-s;\tfrac{\alpha}{\alpha-1})K(s-\gamma).$

This, or slight variants of it, will turn out to have better analytic
properties with respect to $s$.

###  3.3. Check with $\alpha=0$

At $\alpha=0$ the retarded functional equation (3.2) reduces to

(3.8) $H_{0}(s)=\frac{\gamma}{1+\gamma-s}H_{0}(s-\gamma).$

######  Proposition 3.3.

The solution of (3.8) is

(3.9) $H_{0}(s)=\Gamma\left(\frac{1+\gamma-s}{\gamma}\right).$

###### Proof.

Equation (3.8) is a linear, homogeneous first order recurrence relation with
the particular solution (2.3) or (3.9) which also satisfies the boundary
condition $H(1)=1$. It is also the only solution because the only permissible
$\gamma$-periodic function $P(s)$, satisfying $P(s)=P(s-\gamma)$, and making
$H_{0}(s)P(s)$ another solution, is the constant function. This crucial fact
follows because $|H(s)|$ must be bounded (by unity) as $\Im(s)\to\pm\infty$
with $\Re(s)=1$. ∎

We can also verify (3.8) is satisfied by (2.3) from the direct computation of
its Mellin transform as given by (2.2).

###  3.4. Check with $\alpha=1$

For $\alpha=1$ the Gauss hypergeometric function evaluates to

${}_{2}F_{1}(\gamma,1+\gamma-s;2+\gamma-s;1)=\frac{\Gamma(1-\gamma)\Gamma(2+\gamma-s)}{\Gamma(2-s)},$

and the functional equation becomes

(3.10)
$H_{1}(s)=\frac{\gamma\Gamma(1-\gamma)\Gamma(1+\gamma-s)}{\Gamma(2-s)}H_{1}(s-\gamma).$

######  Proposition 3.4 (see Eq. (11), §2.1 of [28]).

The solution of (3.10) is

(3.11)
$H_{1}(s)=\left(\Gamma(1-\gamma)\right)^{(s-1)/\gamma}\frac{\Gamma\left(\frac{\displaystyle
1+\gamma-s}{\displaystyle\gamma}\right)}{\Gamma(2-s)}.$

###### First Proof.

Define $t=s/\gamma$ and $P_{1}(t)$ by

$H_{1}(s)=\left(\Gamma(1-\gamma)\right)^{t}\sin(\pi\gamma t)\Gamma(\gamma
t)\Gamma(\gamma^{-1}-t)P_{1}(t).$

Then (3.10) translates as $P_{1}(t+1)=P_{1}(t)$ and so $P_{1}$ is a
$1$-periodic function. However $P_{1}(t)$ is, in fact a constant, due to
identical arguments given in the proof of Proposition 3.3. This constant value
is

$P_{1}=\frac{1}{\pi\gamma\left(\Gamma(1-\gamma)\right)^{1/\gamma}},$

and the stated result follows. This result can also be found from a
specialisation of Eq. (11), §2.1 of [28]. ∎

###### Second Proof.

We can directly compute the Mellin transform of (2.5) and verify (3.11).
Taking this transform we compute

$H_{1}(s)=\frac{C\,\Gamma(s)}{\pi}\int_{0}^{\infty}dr\;r^{-s}\sin\left(\frac{\pi}{\Gamma(\gamma)}r^{\gamma}\right)e^{-\cos(\pi\gamma)\Gamma(1-\gamma)r^{\gamma}}.$

under the condition $\Re(s)<1+\gamma$ and the more restrictive
$\Re(\gamma)<\frac{1}{2}$, ($\cos(\pi\gamma)>0$). Under an internal change of
variable this integral becomes a known one, [17][pg. 490, 3.944, # 5], with an
evaluation leading to (3.11). In addition we see that $C=1$. ∎

###  3.5. The general solution

Treating the functional equation (3.2) as a recurrence relation and recurring
to the left $k$ times we deduce

(3.12)
$H(s-k\gamma)=\prod_{l=1}^{k}\left(\frac{1+l\gamma-s}{\gamma}\right)\frac{1}{{}_{2}F_{1}(\gamma,1+l\gamma-s;2+l\gamma-s;\alpha)}H(s),$

and thus setting $s=1$

(3.13)
$H_{k}:=H(1-k\gamma)=k!\prod_{l=1}^{k}\frac{1}{{}_{2}F_{1}(\gamma,l\gamma;1+l\gamma;\alpha)}.$

This is meaningful because the denominators never vanish and $H(s)$ is
analytic for $\Re(s)<1+\gamma$ as we demonstrate in the following Lemma.

######  Definition 3.1.

Let $F_{j}:={}_{2}F_{1}(\gamma,j\gamma;1+j\gamma;\alpha)$ for
$j\in\mathbb{Z}_{\geq 0}$.

######  Lemma 3.1.

Let $0<\alpha<1$ and $\gamma>0$. The hypergeometric factors $F_{j}$ satisfy

  1. (1)

$F_{0}=1$,

  2. (2)

are monotonically increasing, $F_{j+1}>F_{j}$, and

  3. (3)

are bounded above, $F_{j}<(1-\alpha)^{-\gamma}$.

###### Proof.

The factors have the integral representation

(3.14)
$F_{j}=(1-\alpha)^{-\gamma}-\alpha\gamma\int_{0}^{1}du\;u^{j\gamma}(1-\alpha
u)^{-\gamma-1}.$

From this we observe that

$\displaystyle
F_{j+1}-F_{j}=\alpha\gamma\int_{0}^{1}du\;u^{j\gamma}(1-u^{\gamma})(1-\alpha
u)^{-\gamma-1}>0.$

In addition $\alpha\gamma\int_{0}^{1}du\;u^{j\gamma}(1-\alpha
u)^{-\gamma-1}>0$ so that $F_{j}<(1-\alpha)^{-\gamma}$. Note that
$F_{\infty}={}_{2}F_{1}(\gamma,-;-;\alpha)=(1-\alpha)^{-\gamma}$. ∎

In contrast recurring to the right encounters the simple poles at
$s=1+k\gamma$, as exemplified by the formula

(3.15)
$H(s+k\gamma)=\gamma^{k}\prod_{l=0}^{k-1}\frac{1}{1-l\gamma-s}{}_{2}F_{1}(\gamma,1-l\gamma-s;2-l\gamma-s;\alpha)\times
H(s).$

######  Proposition 3.5.

The consecutive product of hypergeometric functions has the leading order
behaviour as $N\to\infty$

(3.16)
$\prod_{k=1}^{N}{}_{2}F_{1}(\gamma,1;2+k\gamma-s;\tfrac{\alpha}{\alpha-1})\mathop{\sim}\limits_{N\to\infty}C\left[N+1+\frac{2-s}{\gamma}\right]^{-\alpha/(1-\alpha)},$

where $C$ is a constant independent of $N$. An alternative form of this is

(3.17)
$\prod_{k=1}^{N}{}_{2}F_{1}(\gamma,k\gamma;1+k\gamma;\alpha)\mathop{\sim}\limits_{N\to\infty}C(1-\alpha)^{-N\gamma}\left[N+1+\frac{1}{\gamma}\right]^{-\alpha/(1-\alpha)}.$

###### Proof.

Firstly we write

$(1-\alpha)^{\gamma}{}_{2}F_{1}(\gamma,1-s+k\gamma;2-s+k\gamma;\alpha)={}_{2}F_{1}(\gamma,1;2+k\gamma-s;\tfrac{\alpha}{\alpha-1})=1+A_{k}+B_{k},$

where

$A_{k}=-\frac{\alpha}{1-\alpha}\frac{\gamma}{2-s+k\gamma},\quad
B_{k}=\frac{\alpha^{2}\gamma(\gamma+1)(1-\alpha)^{\gamma}}{2-s+k\gamma}\int_{0}^{1}du\;u^{2-s+k\gamma}(1-\alpha
u)^{-\gamma-2}.$

This is valid because $\Re(s)<1+\gamma$ assures that $2-\Re(s)+k\gamma>0$ for
$k\geq 1$. We will also require $1+A_{k}>0$ and this requires the slightly
stronger condition $2-\Re(s)+k\gamma>\gamma\frac{\alpha}{1-\alpha}$. Splitting
the forgoing sum in the following way

(3.18)
$\log\left(1+A_{k}+B_{k}\right)=\log\left(1+A_{k}\right)+\log\left(1+\frac{B_{k}}{1+A_{k}}\right),$

and employing the elementary inequality
$|\log\left(1+\frac{B_{k}}{1+A_{k}}\right)|\leq\frac{|B_{k}|}{|1+A_{k}|}$
along with the bound

$|B_{k}|\leq\frac{\alpha^{2}}{(1-\alpha)^{2}}\frac{\gamma(\gamma+1)}{(2-s+k\gamma)(3-s+k\gamma)},$

we see that the second term in (3.18) is of ${\rm O}(k^{-2})$ as $k\to\infty$.
Thus in the sum over $k$ it contributes a term of only ${\rm O}(1)$ as
$N\to\infty$.

Now focusing on the first term of (3.18) we see that this is bracketed by
another elementary inequality
$A_{k}-\frac{1}{2}\frac{A_{k}^{2}}{1+A_{k}}\leq\log\left(1+A_{k}\right)\leq
A_{k}-\frac{1}{2}A_{k}^{2}$ (as $A_{k}<0$). Again the second terms of the
lower and upper bounds are of ${\rm O}(k^{-2})$ as $k\to\infty$ and contribute
a term of only ${\rm O}(1)$ as $N\to\infty$. Collecting all of these
contributions in the $k$-sum we have, as $N\to\infty$,

$\displaystyle\sum_{k=1}^{N}\log{}_{2}F_{1}(\gamma,1;2+k\gamma-s;\tfrac{\alpha}{\alpha-1})$
$\displaystyle=-\frac{\alpha}{1-\alpha}\sum_{k=1}^{N}\frac{\gamma}{2-s+k\gamma}+{\rm
O}(1),$ $\displaystyle=-\frac{\alpha}{1-\alpha}\psi(N+1+(2-s)/\gamma)+{\rm
O}(1),$
$\displaystyle=-\frac{\alpha}{1-\alpha}\log\left[N+1+(2-s)/\gamma\right]+{\rm
O}(1).$

Thus (3.16) follows. ∎

For $k\in\mathbb{Z}_{\geq 0}$ we have a denumerable sequence of $K(s_{k})$
values and we know that $K(s)$ is bounded in the finite $s$-plane. Clearly
$K(s)$ is entire, and the question we pose is that of interpolation of an
entire function, i.e. reconstructing it from this sequence of values on
$\mathbb{Z}_{\geq 0}$. This question was addressed in 1884 with the Mittag-
Leffler theorem [37], [10], [7], [46], by the Pringsheim interpolation formula
[44], and by Guichard’s theorem [19]. However use of these involve
restrictions which are not satisfied by $K(s)$.

A more useful function turns out to be

(3.19)
$G(t):=\frac{1}{\Gamma\left(\dfrac{1+\gamma-s}{\gamma}\right)}H(s),\quad
t:=\frac{1-s}{\gamma},$

which eliminates the simple poles at $s=1+k\gamma$, $k\in\mathbb{N}$, so that
in the finite-$t$ plane

(3.20) $|G(t)|\leq(1-\alpha)^{-\gamma}\left[\cosh(\pi\Im(t))\right]^{1/2}.$

However, this is at the expense of a growth $\exp(\pi|\Im(t)|/2)$ as
$\Im(t)\to\pm\infty$. At the prescribed values $s=1-k\gamma$,
$k\in\mathbb{Z}_{\geq 0}$, we deduce

(3.21)
$G_{k}:=G(t=k)=\frac{1}{\prod_{l=1}^{k}{}_{2}F_{1}(\gamma,l\gamma;1+l\gamma;\alpha)},$

which have the growth properties

(3.22)
$G_{k}\mathop{\sim}\limits_{k\to\infty}C(1-\alpha)^{k\gamma}\left[k+1+\frac{1}{\gamma}\right]^{\alpha/(1-\alpha)}.$

######  Definition 3.2.

Recall $F_{j}:={}_{2}F_{1}(\gamma,j\gamma;1+j\gamma;\alpha)$ for
$j\in\mathbb{N}$, and let $t=(1-s)/\gamma$. Define the analytic generating
function $g(x;\gamma,\alpha)$

(3.23)
$g(x):=\sum_{l=0}^{\infty}\frac{(-x)^{l}}{l!}G_{l}=\sum_{l=0}^{\infty}\frac{(-x)^{l}}{l!}\frac{1}{\prod_{j=0}^{l}F_{j}}.$

######  Proposition 3.6.

$g(x)$ is an entire function of $x$ for $\gamma>0$ and $0<\alpha<1$ due to
(3.22).

Because of the bounded growth (3.20) and consequent analyticity of $G(t)$ we
can use the Ramanujan interpolation formula, also known as the Ramanujan
master theorem, given in the first of his quarterly reports of 1913, and
published in [2, §1.2, pg. 298, Ramanujan Notebooks, Part I]

(3.24) $\Gamma(t)G(-t)=\int_{0}^{\infty}dx\;x^{t-1}g(x)=\mathcal{M}[g(x);t].$

See [1, §10.12, pg. 550 and Ex. 34(b)] and [23, Chap. XI] for proofs of this
relation under the same hypotheses of Carlson’s theorem, and [4, App. B] for
recent developments with these methods. Here we recall some details of Hardy’s
theorem and proof, which is framed in terms of the quantity

(3.25) $\phi(u):=\frac{G(u)}{\Gamma(1+u)},$

and the half-plane $\mathcal{H}(\delta):=\\{\Re(u)>-\delta,0<\delta<1\\}$.

######  Proposition 3.7 (Hardy 1940, Chap. XI).

Let

  1. (1)

$\phi(u)$ be regular in $\mathcal{H}(\delta)$,

  2. (2)

$|\phi(u)|<C\exp\left[P\Re(u)+A|\Im(u)|\right]$ in $\mathcal{H}(\delta)$, and

  3. (3)

$A<\pi$.

Then (3.24) is valid.

In our case we have (3.20) and

(3.26)
$|\phi(u)|\leq(1-\alpha)^{-\gamma}\frac{\cosh(\pi\Im(u))}{\Gamma(1+\Re(u))},$

so $A=\pi$, and we need to extend the result of Hardy.

######  Corollary 3.2 ([5], Lemma 3.1).

Prop. 3.7 extends to $A=\pi$, in particular the integrand of the inverse
Mellin transform to (3.24) has algebraic decay

(3.27)
$\left|\Gamma(t)G(-t)\right|\leq(1-\alpha)^{-\gamma}\Gamma(1-\Re(t))\left|\Im(t)\right|^{2\Re(t)-1},$

as $\left|\Im(t)\right|\to\infty$. This means that $G(-t)\in
Q_{\infty}(2,\tfrac{1}{2})$, one of the classes defined in [5]. Consequently
the inverse Mellin transform converges absolutely and uniformly for
$\Re(t)<0$.

###### Proof.

From (3.5) we have

(3.28)
$\left|G(-t)\right|\leq(1-\alpha)^{-\gamma}\frac{\Gamma(1-\Re(t))}{\left|\Gamma(1-t)\right|},$

and from Eq. 2.4.4 of §2.4 of [41] we have

$\displaystyle\log\left|\Gamma(t)\right|\mathop{\to}\limits_{\left|\Im(t)\right|\to\infty}$
$\displaystyle-\tfrac{1}{2}\pi\left|\Im(t)\right|+(\Re(t)-\tfrac{1}{2})\log\left|\Im(t)\right|,$
$\displaystyle\log\left|\Gamma(1-t)\right|\mathop{\to}\limits_{\left|\Im(t)\right|\to\infty}$
$\displaystyle-\tfrac{1}{2}\pi\left|\Im(t)\right|+(\tfrac{1}{2}-\Re(t))\log\left|\Im(t)\right|.$

Thus the exponential decay of the integrand vanishes and one is left with the
algebraic decay of (3.27). ∎

As a consequence of (3.19) we have

(3.29) $\mathcal{M}[h(x);s]=\frac{\Gamma\left(\frac{\displaystyle
1+\gamma-s}{\displaystyle\gamma}\right)}{\Gamma\left(\frac{\displaystyle
s-1}{\displaystyle\gamma}\right)}\mathcal{M}\left[g(x);\frac{s-1}{\gamma}\right].$

This relationship can be translated uniquely into one connecting $h(x)$ and
$g(x)$, given that the inverse Mellin transform of the ratio on the right-hand
side of the above equation is available.

######  Lemma 3.2.

Let $g(x)$ be defined by (3.23). Then the density $h(x;\gamma,\alpha)$ is
given by the Hankel transform of order zero of $g(\cdot)$

(3.30)
$h(x)=\gamma^{2}x^{-\gamma-1}\int_{0}^{\infty}dt\;t^{-\gamma-1}J_{0}(2[xt]^{-\gamma/2})g(t^{-\gamma}).$

Here $J_{0}(z)$ is the standard Bessel function of order zero. The inversion
of this Hankel transform is

(3.31)
$t^{-\gamma/4}g(t^{-\gamma})=\int_{0}^{\infty}dx\;J_{0}(2[tx]^{-\gamma/2})h(x).$

######  Remark 3.1.

A check can be made in the special case of $\alpha=0$. Here $F_{j}=1$ and
$g(x)=e^{-x}$. One finds that (3.30) reduces to

$h(x)=2\gamma
x^{-\gamma-1}\int_{0}^{\infty}dv\;ve^{-v^{2}}J_{0}(2x^{-\gamma/2}v),$

which is a known integral, see [17, pg. 716, 6.631, #(1)], and that this
evaluates to $\frac{1}{2}\exp(-x^{-\gamma})$. Thus we recover $h_{0}(x)$ as
given by (2.2).

######  Remark 3.2.

In the case $\alpha=1$ we find

$F_{j}=\Gamma(1-\gamma)\frac{\Gamma(1+j\gamma)}{\Gamma(1+(j-1)\gamma)},\quad\prod_{j=0}^{l}F_{j}=\left(\Gamma(1-\gamma)\right)^{l}\Gamma(1+l\gamma).$

Consequently we find that the following result.

######  Proposition 3.8.

The generating function $g$ for $\alpha=1$ is a Wright Bessel function of the
form

(3.32) $g(x)=\phi(\gamma,1,-x/\Gamma(1-\gamma)),$

in the notation of [30, pg. 241]. Other notations express this as, see [51],

$g(x)=J^{\gamma}_{0}(x/\Gamma(1-\gamma)),$

or

$g(x)=W(\gamma,1;-x/\Gamma(1-\gamma)).$

Using the known Mellin transform of the Wright Bessel function we can compute
that

$G(t)=\frac{1}{\Gamma\left(\frac{s-1}{\gamma}\right)}\mathcal{M}\left[g(x);\frac{s-1}{\gamma}\right]=\frac{\left(\Gamma(1-\gamma)\right)^{(s-1)/\gamma}}{\Gamma(2-s)}.$

This agrees with the evaluation of $G(t)$ using $H_{1}(s)$ from (3.11).

Combining (3.32) and (3.30) with (2.11) we are lead to the following identity

(3.33)
$\gamma^{2}x^{-\gamma-1}\int_{0}^{\infty}dt\;t^{-\gamma-1}J_{0}(2[xt]^{-\gamma/2})\phi(\gamma,1,-t^{-\gamma}/\Gamma(1-\gamma))\\\
=\gamma^{-1}\left(\Gamma(1-\gamma)\right)^{-1/\gamma}{\rm
H}^{0,1}_{1,1}\left(\frac{x}{\left(\Gamma(1-\gamma)\right)^{1/\gamma}};\begin{array}[]{c}(1-\gamma^{-1},\gamma^{-1})\\\
(0,1)\end{array}\right).$

It remains to directly verify this.

######  Lemma 3.3.

Let $0<\alpha<1$. The finite product $G_{l}$, defined by (3.21), can be
written as an analytic function of complex $\Re(l)>1+\gamma^{-1}$ in an
infinite product form

(3.34) $G_{l}=(1-\alpha)^{\gamma l}\prod_{j=1}^{\infty}\frac{F_{j+l}}{F_{j}},$

and its interpolating function is

(3.35) $G(-t)=(1-\alpha)^{-\gamma
t}\prod_{j=1}^{\infty}\frac{F_{j-t}}{F_{j}}.$

A useful relation, to be used later is

(3.36)
$F_{-t}\prod_{j=1}^{\infty}\frac{F_{j-t}}{F_{j}}=F_{\infty}\prod_{j=0}^{\infty}\frac{F_{j-t}}{F_{j}}.$

###### Proof.

To show convergence of the product, [29, §29, Satz 5, s. 228], we compute
bounds for $\left|\frac{F_{j+l}}{F_{j}}-1\right|$. The factors can be written
(see (3.3))

$F_{j}=(1-\alpha)^{-\gamma}-\alpha\gamma\int_{0}^{1}du\;u^{j\gamma}(1-\alpha
u)^{-\gamma-1}.$

Then we note

$\displaystyle\left|\frac{F_{j+l}}{F_{j}}-1\right|$
$\displaystyle=\frac{\alpha\gamma(1-\alpha)^{\gamma}\int_{0}^{1}du\;(1-\alpha
u)^{-\gamma-1}[u^{j\gamma}-u^{(j+l)\gamma}]}{\left|1-\alpha\gamma(1-\alpha)^{\gamma}\int_{0}^{1}du\;(1-\alpha)^{-\gamma-1}u^{j\gamma}\right|},$
$\displaystyle\leq\frac{\alpha\gamma(1-\alpha)^{-1}\int_{0}^{1}du\;[u^{j\gamma}-u^{(j+l)\gamma}]}{\left|1-\alpha\gamma(1-\alpha)^{-1}\int_{0}^{1}du\;u^{j\gamma}\right|},$
$\displaystyle=\frac{\alpha\gamma^{2}l}{\left[1+(j+l)\gamma\right]\left[1-\alpha(1+\gamma)+j\gamma(1-\alpha)\right]}$
$\displaystyle=\text{${\rm O}(j^{-2})$ as $j\to\infty$},$

uniformly for $l,\alpha,\gamma$ under the conditions stated.

We recall $\lim_{j\to\infty}F_{j}=F_{\infty}=(1-\alpha)^{-\gamma}$ is well
defined. The situation we have is that we take $N\to\infty$ for fixed $l$ so

$\displaystyle\prod_{j=1}^{N}\frac{F_{j+l}}{F_{j}}$
$\displaystyle=\frac{F_{N+1}\ldots F_{N+l}}{F_{1}\ldots F_{l}},\quad\text{as
$N\geq l+1$},$
$\displaystyle\mathop{\to}\limits_{N\to\infty}\frac{1}{F_{1}\ldots
F_{l}}\lim_{N\to\infty}F_{N+1}\ldots F_{N+l},$
$\displaystyle=F_{\infty}^{l}G_{l},$

and (3.34) follows. Now set $l\mapsto l-1$ or $t\mapsto t+1$ in (3.35)

$\displaystyle\prod_{j=1}^{N}\frac{F_{j-1+l}}{F_{j}}$
$\displaystyle=\frac{F_{N+1}\ldots F_{N+l-1}}{F_{1}\ldots
F_{l-1}},\quad\text{as $N\geq l$},$
$\displaystyle\mathop{\to}\limits_{N\to\infty}\frac{1}{F_{1}\ldots
F_{l-1}}\lim_{N\to\infty}F_{N+1}\ldots F_{N+l-1},$
$\displaystyle=F_{\infty}^{l-1}G_{l-1},$
$\displaystyle\prod_{j=1}^{\infty}\frac{F_{j-1+l}}{F_{j}}$
$\displaystyle=\frac{F_{l}}{F_{\infty}}\prod_{j=1}^{\infty}\frac{F_{j+l}}{F_{j}}.$

Thus (3.36) also follows. ∎

In summary we have our penultimate result where the generating function and
density are shown to possess integral representations.

######  Proposition 3.9.

Let $0\leq\alpha<1-\delta$, $\delta>0$. The generating function
$g(x;\gamma,\alpha)$ has the Mellin-Barnes integral representation
$0<c<1+\gamma^{-1}$

(3.37) $g(x)=\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}dt\;x^{-t}\Gamma(t)(1-\alpha)^{-\gamma
t}\prod_{j=1}^{\infty}\frac{{}_{2}F_{1}(\gamma,(j-t)\gamma;1+(j-t)\gamma;\alpha)}{{}_{2}F_{1}(\gamma,j\gamma;1+j\gamma;\alpha)},$

and consequently the density $h(x;\gamma,\alpha)$ has a similar representation
$c<1$

(3.38) $h(x)=\frac{\gamma}{2\pi ix}\int_{c-i\infty}^{c+i\infty}dt\;x^{-\gamma
t}\Gamma(1-t)(1-\alpha)^{-\gamma
t}\prod_{j=1}^{\infty}\frac{{}_{2}F_{1}(\gamma,(j-t)\gamma;1+(j-t)\gamma;\alpha)}{{}_{2}F_{1}(\gamma,j\gamma;1+j\gamma;\alpha)}.$

The Mellin transform of the density $H(s;\gamma,\alpha)$ is

(3.39)
$H(s)=(1-\alpha)^{1-s}\Gamma\left(\frac{1+\gamma-s}{\gamma}\right)\prod_{j=1}^{\infty}\frac{{}_{2}F_{1}(\gamma,1-s+j\gamma;2-s+j\gamma;\alpha)}{{}_{2}F_{1}(\gamma,j\gamma;1+j\gamma;\alpha)}.$

###### Proof.

The exponential generating function sum can be written as a Mellin-Barnes
inversion integral of (3.24)

(3.40) $g(x)=\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}dt\;x^{-t}\Gamma(t)G(-t),$

given convergence of the integrand as $\Im(t)\to\pm\infty$ on the vertical
contour $\Re(t)=c$. The second formula follows from the first under the
absolute convergence of both integrals and the validity of the integral

$\int^{\infty}_{0}du\;u^{1-2t}J_{0}(zu)=2^{1-2t}z^{2t-2}\frac{\Gamma(1-t)}{\Gamma(t)},$

when $\frac{3}{4}<\Re(t)<1$. ∎

###### Second Proof.

There is a simple and direct way to verify (3.38). For generic
$\gamma\notin\mathbb{Q}$ the full set of singularities of the integrand in
(3.38) are: a sequence of simple poles due to $\Gamma(1-t)$ at $t_{k}=k$,
$k\in\mathbb{N}$, and a double sequence of simple poles due to the denominator
factor in ${}_{2}F_{1}(\gamma,(j-t)\gamma;1+(j-t)\gamma;\alpha)$ at
$t_{j,k}=j+\gamma^{-1}k$ for $j,k\in\mathbb{N}$. Thus we are able to displace
the contour to the left so that $c<0$, which will be necessary for our
derivation. Now we substitute (3.38) into the right-hand side of (1.2)

${\rm RHS}=\frac{\gamma^{2}}{2\pi
ix}\int_{c-i\infty}^{c+i\infty}dt\;\Gamma(1-t)(1-\alpha)^{-\gamma
t}\prod_{j=1}^{\infty}\frac{F_{j-t}}{F_{j}}\int_{0}^{x}du\;(x-\alpha
u)^{-\gamma}u^{-1-\gamma t},$

where we have abbreviated the parameter dependence of the hypergeometric
functions and employed the absolute and uniform convergence of the integrals
for $\Re(t)<-1$ to exchange the order. The $u$-integration evaluates to

$\frac{x^{-\gamma-\gamma t}}{-\gamma
t}{}_{2}F_{1}(\gamma,-t\gamma;1-t\gamma;\alpha)=\frac{x^{-\gamma-\gamma
t}}{-\gamma t}F_{-t}.$

Simplifying the right-hand side using (3.36) we find

${\rm RHS}=\frac{\gamma}{2\pi
ix}\int_{c-i\infty}^{c+i\infty}dt\;x^{-\gamma(t+1)}\Gamma(-t)(1-\alpha)^{-\gamma-\gamma
t}\prod_{j=0}^{\infty}\frac{F_{j-t}}{F_{j}},$

and making a translation of the integration variable $v=t+1$ we see

${\rm RHS}=\frac{\gamma}{2\pi
ix}\int_{c+1-i\infty}^{c+1+i\infty}dv\;x^{-\gamma
v}\Gamma(1-v)(1-\alpha)^{-\gamma
v}\prod_{j=0}^{\infty}\frac{F_{j+1-v}}{F_{j+1}}.$

Because $c+1<1$ we can shift the contour one unit to the left thus restoring
the original one, and the product index can be relabelled, we conclude that
the resulting integral is now identical to (3.38). ∎

######  Remark 3.3.

The observations on the large $N$ behaviour of the products in Prop. 3.5 and
the integral representations given above show that the previous Proposition
does not hold in the neighbourhood of $\alpha=1$. When $\alpha=1$ the finite
product has the following evaluation

(3.41)
$\displaystyle\prod_{j=1}^{L}\frac{{}_{2}F_{1}(\gamma,(j-t)\gamma;1+(j-t)\gamma;1)}{{}_{2}F_{1}(\gamma,j\gamma;1+j\gamma;1)}$
$\displaystyle=\frac{1}{\Gamma(1-t\gamma)}\frac{\Gamma(L\gamma+1-t\gamma)}{\Gamma(L\gamma+1)},$
(3.42)
$\displaystyle\mathop{\sim}\limits_{L\to\infty}\frac{1}{\Gamma(1-t\gamma)}\left(L\gamma+\frac{1-t\gamma}{2}\right)^{-t\gamma},$

and so the limit vanishes. Furthermore we know from (3.11) that the inverse
Mellin transform shows, for $c<1+\gamma$,

(3.43) $h_{1}(x)=\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}ds\;x^{-s}\left(\Gamma(1-\gamma)\right)^{(s-1)/\gamma}\frac{\Gamma\left(\frac{\displaystyle
1+\gamma-s}{\displaystyle\gamma}\right)}{\Gamma(2-s)}.$

Clearly there is a transition regime from $\alpha<1$ to $\alpha=1$ where the
limiting behaviour as $L\to\infty$ crosses over from one to the other - a so-
called ”double scaling limit”. We can explicitly demonstrate this with the
following matching procedure: using (3.42) the relevant part of the integrand
of (3.38) is

(3.44) $(1-\alpha)^{-\gamma
t}\prod_{j=1}^{L}\frac{F_{j-t}}{F_{j}}\mathop{\sim}\limits_{L\to\infty\atop\alpha\uparrow
1}\frac{\left[\gamma(1-\alpha)L\right]^{-\gamma t}}{\Gamma(1-t\gamma)}.$

Now if we change variables $t=(s-1)/\gamma$ the integrand of (3.38) will match
that of (3.43) if we simultaneously take $L\to\infty$ with $\alpha\to
1^{\uparrow}$ such that the product

$(1-\alpha)L=\frac{1}{\gamma\left(\Gamma(1-\gamma)\right)^{1/\gamma}},$

is fixed.

###  3.6. The large $\gamma$ case

A particularly instructive case of our foregoing analysis is when $\gamma$ is
real and large whilst $0<\alpha<1$. In this way we have a simple example where
$\alpha$ is explicitly exhibited. It is apparent that in the opposite case of
$\gamma\to 0$ only the trivial solution of $h(x)$ with Dirac-atomic measure at
the origin is admissible. In this regime the leading order behaviour of the
hypergeometric function is [40], [9] is

(3.45) ${}_{2}F_{1}(a+\epsilon\lambda,b;c+\lambda;z)\sim\left(1-\epsilon
z\right)^{-b},\quad|\lambda|\to\infty,$

excluding the region around the singular point $z_{S}=1/\epsilon$. In our
application we require the connection formula

(3.46)
${}_{2}F_{1}(\gamma,(k-t)\gamma;1+(k-t)\gamma;\alpha)=(1-\alpha)^{-\gamma}{}_{2}F_{1}(\gamma,1;1+(k-t)\gamma;-\tfrac{\alpha}{1-\alpha}),$

with $\lambda=(k-t)\gamma$ and $0<\epsilon=1/(k-t)<1$ so our singular point is

(3.47) $\alpha_{S}=1+1/(k-1-t)>1,$

which is excluded from our consideration. Using the definition 3.1 we have

(3.48)
$F_{j}\mathop{\sim}\limits_{\gamma\to\infty}(1-\alpha)^{-\gamma}\left[1+\frac{1}{j}\frac{\alpha}{1-\alpha}\right]^{-1},\quad
j\geq 1.$

With this knowledge one can easy verify the following results.

######  Proposition 3.10.

For $0\leq\alpha<1$ and for $\gamma\to\infty$ the leading order hypergeometric
product $G_{l}$, $l\geq 0$ is

(3.49)
$G_{l}=(1-\alpha)^{l\gamma}\frac{\Gamma(l+\frac{1}{1-\alpha})}{l!\Gamma(\frac{1}{1-\alpha})}.$

Consequently the generating function is given by the confluent hypergeometric
function

(3.50) $g(x)={}_{1}F_{1}(\tfrac{1}{1-\alpha};1;-(1-\alpha)^{\gamma}x).$

This is an entire function of $x$ and $\frac{1}{1-\alpha}$ with exponential
decay as $x\to\infty$ and its Mellin transform is

(3.51) $\mathcal{M}[g(x);t]=(1-\alpha)^{-\gamma
t}\frac{\Gamma(t)\Gamma(\frac{1}{1-\alpha}-t)}{\Gamma(1-t)\Gamma(\frac{1}{1-\alpha})},$

for $0<\Re(t)<\frac{1}{1-\alpha}$, using Oberhettinger’s tables [38, 13.49,
pg. 145]. Using the Master Theorem (3.24) we obtain the interpolating function

(3.52) $G(t)=(1-\alpha)^{\gamma
t}\frac{\Gamma(\frac{1}{1-\alpha}+t)}{\Gamma(1+t)\Gamma(\frac{1}{1-\alpha})}.$

###### Proof.

The steps indicated above are self-explanatory. One can also directly solve
the functional equation (3.2) directly. Expanding to sufficient order for
large $\gamma$

(3.53)
$\frac{\gamma}{1+\gamma-s}{}_{2}F_{1}(\gamma,1+\gamma-s;2+\gamma-s;\alpha)\\\
=\frac{\gamma}{1+\gamma-s}(1-\alpha)^{1-\gamma}{}_{2}F_{1}(2-s,1;2+\gamma-s;\alpha),\\\
=(1-\alpha)^{1-\gamma}\left\\{1+\frac{\alpha}{\gamma}+\frac{(\alpha-1)(1-s)}{\gamma}+\rm{O}(\gamma^{-2})\right\\},\\\
=(1-\alpha)^{-\gamma}\left[\frac{1-\alpha/\gamma}{1-\alpha}+\frac{1-s}{\gamma}\right]^{-1}\left\\{1+\rm{O}(\gamma^{-2})\right\\}.$

Solving the functional equation

(3.54)
$\frac{H(s)}{H(s-\gamma)}=\frac{(1-\alpha)^{-\gamma}}{\frac{1-\alpha/\gamma}{1-\alpha}+\frac{1-s}{\gamma}},$

in the same manner as in the proof of (3.4) and appealing to the bound (3.5)
for $\Im(s)\to\pm\infty$ we deduce

(3.55)
$H(s)=(1-\alpha)^{1-s}\frac{\Gamma(\frac{1-\alpha/\gamma}{1-\alpha}+\frac{1-s}{\gamma})}{\Gamma(\frac{1-\alpha/\gamma}{1-\alpha})}.$

Recalling (3.29) this gives us (3.52) after discarding the small term
$\alpha/\gamma$. ∎

######  Remark 3.4.

Inverting the Mellin transform (3.55) gives an elementary function form for
the density

(3.56)
$h_{\gamma}(x)=\frac{(1-\alpha)^{-\frac{\gamma-\alpha}{1-\alpha}}}{\Gamma\left(\frac{1-\alpha/\gamma}{1-\alpha}\right)}x^{-1-\frac{\gamma-\alpha}{1-\alpha}}\exp\left[-(1-\alpha)^{-\gamma}x^{-\gamma}\right].$

###  3.7. Continued Fraction Forms

Our final representation of $g(x)$ is as a continued fraction, however this is
best expressed in terms of a related function rather than $g$ itself.

######  Definition 3.3.

Let $L(z;\gamma,\alpha)$ be the Laplace transform of the generating function
$g(x)$, for $\Re(z)>0$.

(3.57)
$L(z):=\mathcal{L}\left[g(x);z\right]=\int_{0}^{\infty}dx\;e^{-zx}g(x).$

Thus we have the convergent expansion due to (3.22) for
$|z|>(1-\alpha)^{\gamma}$

(3.58)
$L(z)=\sum_{l=0}^{\infty}\frac{(-1)^{l}}{\prod_{j=0}^{l}F_{j}}z^{-l-1}.$

######  Proposition 3.11.

$L(z)$ is a special $M$-continued fraction, see [8, pg. 38, Chap. 1], given by

(3.59)
$\frac{1}{L(z)}=z+\vbox{\hbox{\huge$\mathrm{K}$}}_{m=1}^{\infty}\cfrac{\displaystyle
F_{m-1}z}{\displaystyle F_{m}z-1},$

where our notation is equivalent to the alternative forms

$b_{0}+\vbox{\hbox{\huge$\mathrm{K}$}}^{\infty}_{j=1}\frac{a_{j}}{b_{j}}=b_{0}+\cfrac{a_{1}}{b_{1}+\cfrac{a_{2}}{b_{2}+\cfrac{a_{3}}{b_{3}+\ldots}}}=b_{0}+\frac{a_{1}}{b_{1}}\genfrac{}{}{0.0pt}{}{}{+}\frac{a_{2}}{b_{2}}\genfrac{}{}{0.0pt}{}{}{+}\frac{a_{3}}{b_{3}}\genfrac{}{}{0.0pt}{}{}{+\ldots}.$

###### Proof.

Eq. (3.59) is found by adapting a result known to Euler, see [8, (1.7.2)], for
the recasting of a series

(3.60) $e_{n}=1+\sum_{k=1}^{n}\prod_{j=1}^{k}\big{(}-E_{j}\big{)}z^{k},$

into

$\displaystyle
e_{n}=\cfrac{1}{1+\vbox{\hbox{\huge$\mathrm{K}$}}^{n}_{j=1}\cfrac{\displaystyle
E_{j}z}{\displaystyle 1-E_{j}z}}.$

It is straight forward to show the partial numerators and denominators $A_{n}$
and $B_{n}$ of $e_{n}=\frac{A_{n}}{B_{n}}$ satisfy the recurrence relations

$\displaystyle A_{n}$ $\displaystyle=b_{n}A_{n-1}+a_{n}A_{n-2},$
$\displaystyle B_{n}$ $\displaystyle=b_{n}B_{n-1}+a_{n}B_{n-2},$

with $b_{0}=0,a_{1}=1,b_{1}=1$ and for $m\geq 1$, $a_{m+1}=E_{m}z$ and
$b_{m+1}=1-E_{m}z$. The solutions to this recurrence give $B_{n}=1$ and
$A_{n}$ the right-hand side of (3.60). ∎

######  Remark 3.5.

As a check for $\alpha=0$ we recall $F_{j}=1$, $j\geq 0$ so that

$\frac{1}{L_{0}(z)}=z+\cfrac{z}{z-1+\cfrac{z}{z-1+\cfrac{z}{\ddots}}}=z+R,$

however we observe $R=\dfrac{z}{z-1+R}$ so that $(R-1)(R+z)=0$. Choosing $R=1$
as the appropriate root we verify that

(3.61) $L_{0}(z)=\frac{1}{z+1},$

which should be compared with the Laplace transform of $g_{0}(x)=e^{-x}$.

######  Remark 3.6.

For $\alpha=1$ we can record a result which may be known

(3.62)
$\frac{1}{L_{1}(z)}=z+\vbox{\hbox{\huge$\mathrm{K}$}}_{m=1}^{\infty}\cfrac{\frac{\Gamma(1+(m-1)\gamma)}{\Gamma(1+(m-2)\gamma)}\Gamma(1-\gamma)z}{\frac{\Gamma(1+m\gamma)}{\Gamma(1+(m-1)\gamma)}\Gamma(1-\gamma)z-1}.$

On the other hand the Laplace transform of $g_{1}$ is the Mittag-Leffler
function $E_{\gamma}(\cdot)$, see [30, pg. 269]

(3.63)
$L_{1}(z)=\frac{1}{z}E_{\gamma}\big{(}-\tfrac{1}{\Gamma(1-\gamma)z}\big{)}.$

This is an entire function of $z$ of order $\frac{1}{\Re(\gamma)}$ and type
unity.

######  Remark 3.7.

From the results of Prop. 3.10 we can show that for $\gamma\to\infty$ the
Laplace transform has the algebraic form

(3.64)
$L_{\gamma}(z)=z^{-1}\left(1+(1-\alpha)^{\gamma}z^{-1}\right)^{-\frac{1}{1-\alpha}}.$

This can be established from the convergent expansion (3.58) about $z=\infty$
using (3.49) or computing the Laplace transform of (3.50) with the assistance
of [12, pg. 215, §2.2 (11)] assuming $\Re(z)>0$. Therefore we note the
M-continued fraction form

(3.65)
$\frac{1}{L_{\gamma}(z)}=z+\vbox{\hbox{\huge$\mathrm{K}$}}_{m=1}^{\infty}\cfrac{\frac{(m-1)}{(m-1+\frac{\alpha}{1-\alpha})}(1-\alpha)^{-\gamma}z}{\frac{m}{(m+\frac{\alpha}{1-\alpha})}(1-\alpha)^{-\gamma}z-1}.$

##  4\. Numerical Evaluations

The most robust numerical strategy for computing the density is to use the
Mellin-Barnes integral representation (3.38) because the integrand is smooth
and decays exponentially as $\Im(t)\to\pm\infty$. There are oscillations but
they have a period of the same order as the foregoing decay length and
therefore are rapidly damped. Our methods employed 64-bit or real (kind=8)
precision with the gfortran F90 compiler GNU Fortran (Ubuntu 7.5.0-3ubuntu1
18.04) 7.5.0 and software from the following sources:

  * •

We have employed DQAGI for the Mellin-Barnes integrals and DQAGS for the
normalisation and integral equations from the QUADPACK package[43]. The
absolute and relative error goals were $10^{-6}$ and $10^{-4}$ respectively. A
convenient choice for $c$ is half.

  * •

The complex gamma function was computed with the error control ACM algorithm
421 [31] with the error value set at $10^{-8}$.

  * •

The complex Gauss hypergeometric function was computed using algorithm
AEAE_v1_0 from the CPC program library [36]. An error criteria is computed
from certain contiguous relations.

  * •

Extrapolation of the Gaussian hypergeometric function products up to $2^{13}$
factors was performed by ACM algorithm 602 HURRY, an acceleration algorithm
for scalar sequences and series [13].

Two criteria were used to quantify the all-up global errors - the
normalisation of the density and the residual of the solution in satisfying
the integral equation (1.2) at $x=1$. These are tabulated in Tables 1 and 2
respectively over a range of $\gamma\in(0,1]$ and $\alpha\in(0,1)$ values. Our
error settings are relatively weak and no attempt has been made in pushing the
methods to their accuracy limits.

Table 1. Absolute Normalisation Errors $\gamma$ $\alpha$ | 0.10 | 0.20 | 0.30 | 0.40 | 0.50 | 0.60
---|---|---|---|---|---|---
0.25 | 0.318E-05 | 0.254E-04 | 0.633E-04 | 0.333E-06 | 0.888E-05 | 0.111E-05
0.50 | 0.331E-05 | 0.126E-05 | 0.201E-06 | 0.122E-04 | 0.733E-05 | 0.138E-05
0.75 | 0.162E-04 | 0.516E-05 | 0.301E-05 | 0.128E-05 | 0.106E-06 | 0.150E-05
1.00 | 0.160E-04 | 0.389E-04 | 0.664E-05 | 0.380E-06 | 0.260E-05 | 0.157E-04
Table 2. Absolute Residual Errors $\gamma$ $\alpha$ | 0.10 | 0.20 | 0.30 | 0.40 | 0.50 | 0.60
---|---|---|---|---|---|---
0.25 | 0.237E-06 | 0.143E-07 | 0.653E-07 | 0.209E-06 | 0.689E-09 | 0.746E-07
0.50 | 0.126E-05 | 0.517E-08 | 0.122E-06 | 0.176E-06 | 0.677E-07 | 0.341E-07
0.75 | 0.622E-05 | 0.347E-05 | 0.284E-05 | 0.383E-06 | 0.253E-06 | 0.129E-06
1.00 | 0.596E-06 | 0.369E-05 | 0.131E-06 | 0.485E-06 | 0.301E-05 | 0.876E-05

We present our results in a sequence of Figures 2, 3, 4, 5, 5 and 6 plotting
the density $h(x;\gamma,\alpha)$ versus $x$ with increasing values of $\gamma$
and within each plot $\alpha$ increases through the range $[0,1)$.

Figure 2. Densities $h(x;\tfrac{1}{4},\alpha)$ versus $x$ for
$\alpha=\\{0,\frac{1}{4},\frac{1}{2},\frac{3}{4},1\\}$ computed using (3.38)
and (2.5). Figure 3. Densities $h(x;\tfrac{1}{2},\alpha)$ versus $x$ for
$\alpha=\\{0,\frac{1}{4},\frac{1}{2},\frac{3}{4},1\\}$ computed using (3.38)
and (2.5).  Figure 4. Densities $h(x;\tfrac{3}{4},\alpha)$ versus $x$ for
$\alpha=\\{0,\frac{1}{4},\frac{1}{2},\frac{3}{4},1\\}$ computed using (3.38)
and (2.5).  Figure 5. Densities $h(x;1,\alpha)$ versus $x$ for
$\alpha=\\{0,\frac{1}{4},\frac{1}{2},\frac{3}{4}\\}$ computed using (3.38).

Figure 6. Large-$\gamma$ densities as computed using (3.56). In this case
$\gamma=4$ and
$\alpha\in\\{0,\tfrac{1}{4},\tfrac{1}{2},\tfrac{3}{4},\tfrac{7}{8}\\}$.

The densities are unimodal and exhibit a tail as $x\to\infty$ with algebraic
decay, having the leading form $x^{-1-\gamma}$ for all $\gamma>0$ and
$0<\alpha\leq 1$. For non-rational $\gamma$ there is a cascade of lower order
terms with algebraic decay $x^{-1-j\gamma}$ for $j\in\mathbb{N}$. This is
illustrated by formula (2.7) in the case $\alpha=1$ and from the general
formula (3.38) for $0<\alpha<1$. This latter formula can be evaluated as a sum
over the algebraically decaying terms by enclosing the poles of the Gamma
function at $t\in\mathbb{N}$ to the right of the contour, similar to (2.7)

(4.1) $h(x)=\gamma
x^{-1}\sum_{j=0}^{\infty}\frac{(-)^{j}}{j!}\left[(1-\alpha)x\right]^{-\gamma(j+1)}\prod_{k=0}^{j}{}_{2}F_{1}(\gamma,-k\gamma;1-k\gamma;\alpha).$

For fixed $\gamma$ and $\alpha$ increasing from $0$ to $1$ the maxima occurs
at increasing $x$-values and the maximal density itself decreases. The
locations and densities of the maxima at $\alpha=0$ are easily seen to be

(4.2) $x_{max}=\left(\frac{\gamma+1}{\gamma}\right)^{-1/\gamma},\quad
h_{max}=e^{-\frac{1}{\gamma}-1}\left(\frac{1}{\gamma}+1\right)^{1/\gamma}(\gamma+1).$

This latter situation is where the bounds (3.4) are sharp.

Table 3. Maxima for $\alpha=1$ and some rational $\gamma$ $\gamma$ | $x_{max}$ | $h_{max}$
---|---|---
1/4 | 0.00539575 | 2.61146
1/3 | 0.0579035 | 0.744108
1/2 | 0.523599 | 0.294463
2/3 | 1.70192 | 0.217484
3/4 | 2.82353 | 0.202068

The behaviour as $x\to 0^{+},{\rm arg}(x)=0$ is quite different, however. The
density has an algebraic divergence here, but it is totally overwhelmed by an
exponential suppression. This is simply illustrated by the large-$\gamma$
formula (3.56) and the $\alpha=1$ special cases (2.8). The asymptotic
behaviour of the Fox H-function as $x\to 0$ with ${\rm arg}(x)=0$ is known,
see [35] pg. 2 Definition 1.1 and Theorem 1.1, Case 4, Eq. (1.109) in (ii) of
Theorem 1.3, so for $0<\gamma<1$ and $\alpha=1$ we have

(4.3) $h_{1}(x)\mathop{\sim}\limits_{x\to
0^{+}}\frac{1}{\gamma\left(\Gamma(1-\gamma)\right)^{1/\gamma}}x^{-(1-\frac{1}{2}\gamma)/(1-\gamma)}\\\
\times\exp\left[-(1-1/\gamma)\gamma^{1/(1-\gamma)}\left(\Gamma(1-\gamma)\right)^{1/(1-\gamma)}x^{-\gamma/(1-\gamma)}\right].$

An outstanding task is to find the analogous leading order asymptotics for
$x\to 0$ and general $\alpha$.

##  5\. Discussion

We pick up here the thread in Section §1 concerning related works in the
probability literature. Schlather [48] treats a different operation from the
one considered in [18], [25] or here in that he takes a $1/p$ root of a sum of
$p^{\rm th}$ powers. In addition to being an altogether different problem
producing a very different solution, he employs a different method. For
example he finds a different family of limit distributions parametrised by $c$
with $0\leq c\leq\infty$. And he discusses domains of attraction, but with
norming depending on $c$ (unlike in the foregoing works) and of course
different limits.

One property of $h(x;\gamma,\alpha)$ for fixed $\alpha$ and $\gamma$ is that
because it is obtained as a normed limit it is infinitely divisible, and so
has a Lévy measure, in the Kolmogorov or Lévy canonical representations of the
characteristic function. Can we see this transparently from our product
representation for $H(s)$, as given by (3.39)? The distribution
$h(x;\gamma,\alpha)$ is also associated with a Lévy process, i.e. a process
with stationary, independent increments. For $\alpha=1$ we have recovered a
known identification of $h(x;\gamma,\alpha)$ as a Fox $H$-function, and we
have mentioned some previous works [49], [39], [28] where there is a nice
statistical interpretation of this fact. However the sequence (1.1) is in the
domain of attraction of the normal distribution for $\alpha=1$, so there is a
”discontinuity”: (1.1) interpolates between Fréchet and normal, whereas (1.2)
interpolates between Fréchet and this Fox $H$-function.

Regarding the nature of the generic solution $h(x;\gamma,\alpha)$ we have
found here, it appears that this is a novel extension of the hypergeometric
function to the best of our knowledge. Standard extensions through the Mittag-
Leffler, Wright-Bessel and Meijer-$G$ functions to the Fox H-functions clearly
appear when $\alpha=1$, but not for $0<\alpha<1$. The clearest way to
distinguish the generic functions found here is to think of the standard
hypergeometric or Fox-H function as series

$H\left(\begin{array}[]{ccc}(A_{1},a_{1})&\ldots&(A_{P},a_{P})\\\
(B_{1},b_{1})&\ldots&(B_{Q},b_{Q})\end{array};x\right)=\sum_{n\geq
0}C_{n}x^{n},$

and that the coefficients $C_{n}$ satisfy a first order difference equation
with respect to $n$, in unit steps in the forward direction $n\mapsto n+1$,

$C_{n+1}=\frac{(A_{1}+a_{1}n)\ldots(A_{P}+a_{P}n)}{(B_{1}+b_{1}n)\ldots(B_{Q}+b_{Q}n)}C_{n},\left\\{A_{1},a_{1},\ldots,A_{P},a_{P},B_{1},b_{1},\ldots,B_{Q},b_{Q}\right\\}\in\mathbb{C},$

i.e. rational coefficients in $n$. Such equations can be solved in terms of
products of gamma functions or Pochhammer symbols given suitable boundary
conditions. However in our generic case the coefficient of the first order
difference equation is itself a hypergeometric function of $n$, which is the
solution of a second-order difference equation with respect to $n$, but as
$n\mapsto n+1/\gamma$. Alternatively the original coefficients $C_{n}$, or
effectively $H(s=1-n\gamma)$, would satisfy a non-linear difference equation
with rational coefficients in $n$ but with incommensurate steps of $1$ and
$\gamma^{-1}$

(5.4)
$0=\left[1+\gamma-s+\alpha(1-s)\right]\frac{H(s)}{H(s-\gamma)}+\alpha(s-2)\frac{H(s-1)}{H(s-1-\gamma)}+(s-\gamma)\frac{H(s+1)}{H(s+1-\gamma)}.$

This is an entirely higher level of special function.

It may be of interest to note that another family of extensions to the Gauss
hypergeometric functions has been explored recently, under the monikers of
extended hypergeometric functions [6], extended $\tau$-Gauss hypergeometric
functions [50], [42], extended, generalised $\tau$-Gauss hypergeometric
functions [32]. The common theme here is that a pair of Pochhammer symbol
ratios in $C_{n}$, which are equivalent to the beta function $B(a,b)$ [45,
§5.12], are replaced by

(5.5)
$B(a,b;\sigma)=\int_{0}^{1}dt\;t^{a-1}(1-t)^{b-1}\exp\left[-\frac{\sigma}{t(1-t)}\right],\Re(\sigma)>0,$

along with

$F_{\sigma}(a,b;c;z)=\sum_{n=0}^{\infty}(a)_{n}\frac{B(b+n,c-b;\sigma)}{B(b,c-b)}\frac{z^{n}}{n!},$

or ultimately

$B(a,b;\sigma;\alpha,\beta,\rho,\lambda)=\int_{0}^{1}dt\;t^{a-1}(1-t)^{b-1}{}_{1}F_{1}\left(\alpha;\beta;-\frac{\sigma}{t^{\rho}(1-t)^{\lambda}}\right),\Re(\sigma)>0.$

There is no evidence that the extended beta function (5.5) satisfies a
difference equation like that of (5.4).

##  6\. Acknowledgements

The authors would like to acknowledge Jesse Goodman for organising the 2020
Taupo Probability Workshop and the funding support of the University of
Auckland, which facilitated this collaboration.

## References

  * [1] G. E. Andrews, R. Askey, and R. Roy.  Special functions, volume 71 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 1999.
  * [2] B. C. Berndt.  Ramanujan’s notebooks. Part I.  Springer-Verlag, New York, 1985.  With a foreword by S. Chandrasekhar.
  * [3] L. Bogachev, G. Derfel, S. Molchanov, and J. Ockendon.  On bounded solutions of the balanced generalized pantograph equation.  In Topics in stochastic analysis and nonparametric estimation, volume 145 of IMA Vol. Math. Appl., pages 29–49. Springer, New York, 2008\.
  * [4] B. Candelpergher.  Ramanujan summation of divergent series, volume 2185 of Lecture Notes in Mathematics.  Springer, Cham, 2017.
  * [5] M. A. Chaudhry and A. Qadir.  Extension of Hardy’s class for Ramanujan’s interpolation formula and master theorem with applications.  J. Inequal. Appl., pages 2012:52, 13, 2012.
  * [6] M. Aslam Chaudhry, A. Qadir, M. Rafique, and S. M. Zubair.  Extension of Euler’s beta function.  J. Comput. Appl. Math., 78(1):19–32, 1997.
  * [7] J. B. Conway.  Functions of one complex variable, volume 11 of Graduate Texts in Mathematics.  Springer-Verlag, New York-Berlin, second edition, 1978.
  * [8] V. B. Cuyt, A.and Petersen, B. Verdonk, H. Waadeland, and W. B. Jones.  Handbook of continued fractions for special functions.  Springer, New York, 2010.
  * [9] M. Cvitković, A.-S. Smith, and J. Pande.  Asymptotic expansions of the hypergeometric function with two large parameters—application to the partition function of a lattice gas in a field of traps.  J. Phys. A, 50(26):265206, 24, 2017.
  * [10] P. J. Davis.  Interpolation and approximation.  Dover Publications, Inc., New York, 1975.  Republication, with minor corrections, of the 1963 original, with a new preface and bibliography.
  * [11] G. A. Derfel.  A probabilistic method for studying a class of functional-differential equations.  Ukrain. Mat. Zh., 41(10):1322–1327, 1436, 1989.
  * [12] A. Erdélyi, W. Magnus, F. Oberhettinger, and F. G. Tricomi.  Tables of integral transforms. Vol. I.  McGraw-Hill Book Company, Inc., New York-Toronto-London, 1954.
  * [13] T. Fessler, D. A. Smith, and W. Ford.  Algorithm 602: Hurry: An acceleration algorithm for scalar sequences and series.  ACM Trans. Math. Softw., 9(3):355–357, September 1983.
  * [14] R.A. Fisher and L.H.C. Tippett.  Limiting forms of the frequency distribution of the largest and smallest member of a sample.  Proc. Cambridge Philosophical Society, 24:180–190, 1928.
  * [15] M. Fréchet.  Sur la loi de probabilité de l’écart maximum.  Ann. Soc. Polon. Math., 6:93, 1927.
  * [16] D. P. Gaver, Jr.  An absorption probability problem.  J. Math. Anal. Appl., 9:384–393, 1964.
  * [17] I. S. Gradshteyn and I. M. Ryzhik.  Table of integrals, series, and products.  Academic Press [Harcourt Brace Jovanovich, Publishers], New York-London-Toronto, Ont., 1980.
  * [18] P. E. Greenwood and G. Hooghiemstra.  On the domain of attraction of an operator between supremum and sum.  Probab. Theory Related Fields, 89(2):201–210, 1991.
  * [19] C. Guichard.  Sur les fonctions entières.  Annales scientifiques de l’École Normale Supérieure, 3e série, 1:427–432, 1884.
  * [20] E. J. Gumbel.  Statistics of extremes.  Columbia University Press, New York, 1958.
  * [21] A. J. Hall and G. C. Wake.  A functional-differential equation arising in modelling of cell growth.  J. Austral. Math. Soc. Ser. B, 30(4):424–435, 1989.
  * [22] A. J. Hall and G. C. Wake.  Functional differential equations determining steady size distributions for populations of cells growing exponentially.  J. Austral. Math. Soc. Ser. B. Appl. Math., 31(4):434–453, 1990\.
  * [23] G. H. Hardy.  Ramanujan. Twelve lectures on subjects suggested by his life and work.  Cambridge University Press, Cambridge, England; Macmillan Company, New York, 1940.
  * [24] P. Hernández-Llanos, Y. Quintana, and A. Urieles.  About extensions of generalized apostol-type polynomials.  Results in Mathematics, 68(1):203–225, 2015.
  * [25] G. Hooghiemstra and P. E. Greenwood.  The domain of attraction of the $\alpha$-sun [$\alpha$-sum] operator for type II and type III distributions.  Bernoulli, 3(4):479–489, 1997.
  * [26] A. Iserles and Y. Liu.  Integro-differential equations and generalized hypergeometric functions.  J. Math. Anal. Appl., 208(2):404–424, 1997.
  * [27] A. Janßen.  On Some Connections between Light Tails, Regular Variation and Extremes.  PhD thesis, 2010.
  * [28] W. Jedidi, T. Simon, and M. Wang.  Density solutions to a class of integro-differential equations.  J. Math. Anal. Appl., 458(1):134–152, 2018.
  * [29] K. Knopp.  Theorie and Anwendung der unendlichen Reihen.  Fünfte berichtigte Auflage. Die Grundlehren der Mathematischen Wissenschaften, Band 2. Springer-Verlag, Berlin-New York, 1964.
  * [30] A. Kochubei and Y. Luchko, editors.  Handbook of Fractional Calculus with Applications, volume 1. Basic Theory of De Gruyter Reference.  De Gruyter, 2019.
  * [31] H. Kuki.  Complex gamma function with error control.  Commun. ACM, 15(4):262–267, April 1972.
  * [32] D. Kumar, R. K. Gupta, and B. S. Shaktawat.  Certain results on extended generalized $\tau$-Gauss hypergeometric function.  Honam Math. J., 38(4):739–752, 2016.
  * [33] J. C. Lagarias and W.-C. W. Li.  The Lerch zeta function I. Zeta integrals.  Forum Mathematicum, 24(1):1, 2012.
  * [34] Q.-M. Luo.  Apostol-Euler polynomials of higher order and Gaussian hypergeometric functions.  Taiwanese J. Math., 10(4):917–925, 06 2006.
  * [35] A. M. Mathai, R. K. Saxena, and H. J. Haubold.  The H-Function : Theory and Applications.  New York : Springer, 2010.
  * [36] N. Michel and M. V. Stoitsov.  Fast computation of the Gauss hypergeometric function with all its parameters complex with application to the Pöschl-Teller-Ginocchio potential wave functions.  Comput. Phys. Comm., 178(7):535–551, 2008.
  * [37] G. Mittag-Leffler.  Sur la représentation analytique des fonctions monogènes uniformes.  Acta Math., 4(1):1–79, 1884.  D’une variable indépendante.
  * [38] F. Oberhettinger.  Tables of Mellin transforms.  Springer-Verlag, New York, 1974.
  * [39] A. G. Pakes.  On generalized stable and related laws.  J. Math. Anal. Appl., 411(1):201–222, 2014.
  * [40] R. B. Paris.  Asymptotics of the Gauss hypergeometric function with large parameters, I.  J. Class. Anal., 2(2):183–203, 2013.
  * [41] R. B. Paris and D. Kaminski.  Asymptotics and Mellin-Barnes integrals, volume 85 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 2001.
  * [42] R. K. Parmar.  Extended $\tau$-hypergeometric functions and associated properties.  C. R. Math. Acad. Sci. Paris, 353(5):421–426, 2015.
  * [43] R. Piessens, E. de Doncker-Kapenga, C. W. Überhuber, and D. K. Kahaner.  QUADPACK, volume 1 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 1983.
  * [44] A. Pringsheim.  Vorlesungen über funktionenlehre: Eindeutige analytische funktionen. II, 2.  Mathematischen Wissenschaften mit einschluss ihrer anwendungen. Verlag und Druck von B. G. Teubner, Leipzig, 1932.
  * [45] NIST Digital Library of Mathematical Functions.  http://dlmf.nist.gov/, Release 1.0.9 of 2014-08-29.
  * [46] R. Remmert.  Classical topics in complex function theory, volume 172 of Graduate Texts in Mathematics.  Springer-Verlag, New York, 1998.  Translated from the German by Leslie Kay.
  * [47] S. I. Resnick.  Extreme values, regular variation, and point processes, volume 4 of Applied Probability. A Series of the Applied Probability Trust.  Springer-Verlag, New York, 1987.
  * [48] M. Schlather.  Limit distributions of norms of vectors of positive i.i.d. random variables.  Ann. Probab., 29(2):862–881, 2001.
  * [49] W. R. Schneider.  Generalized one-sided stable distributions.  In Stochastic processes—mathematics and physics, II (Bielefeld, 1985), volume 1250 of Lecture Notes in Math., pages 269–287. Springer, Berlin, 1987.
  * [50] N. Virchenko, S. L. Kalla, and A. Al-Zamel.  Some results on a generalized hypergeometric function.  Integral Transform. Spec. Funct., 12(1):89–100, 2001.
  * [51] E. M. Wright.  The asymptotic expansion of the generalized hypergeometric function.  J. London Math. Soc., 10:287–293, 1935.