\section{Introduction}

Confinement of light to wavelength-scale volumes possesses important
implications to control the light-matter interactions at single photon and
single emitter level. The emerging field of nanofabrication technology has
encouraged various designs of micro/nano-photonic structures to be
investigated. In particular 2-D \cite{scherer pcm} and 1-D \cite{foresi,
painter nb} photonic crystal (PhC) structures have enabled new prospects in
nanolasers \cite{painter}, chemical sensing \cite{loncar sensing} and
optical switching \cite{tanabe}. Moreover, the strong light-matter coupling
\cite{nature pcm, imamoglu pcm, dirk, lukin pcm} in these PhC structures,
has opened new avenues for quantum information technology. However,
integrating these PhC structures to fiber-based communication channels
still remains challenging. In this context, tapered optical fibers with subwavelength diamters, known
as optical nanofibers, may provide a promising platform for realizing strong light-matter
interaction in a fiber-based system.

In recent years, optical nanofibers have
opened a new paradigm in manipulating single/few-atom fluorescence \cite%
{opex, njp, chandra opex}. A nanofiber is realized by adiabatically tapering
a conventional single mode optical fiber which enables light in the
fundamental mode to be coupled in and out of the nanofiber with more than $%
95\%$ efficiency. Due to the subwavelength diameter the field in the guided
mode is strongly confined in the transverse direction. As a result the
spontaneous emission of atoms can be strongly modified around the nanofiber
and an appreciable amount of atomic fluorescence can be coupled to the
guided modes \cite{spont em, chandra prl}. Also the strong evanescent tail
of the guided field gives excellent access for trapping and probing few
atoms around the nanofiber \cite{scattering, trap1, arno, kimble trap}.

The light-matter coupling in such nanofiber system can be substantially
improved by fabricating a 1-D PhC structure on the nanofiber. The PhC
structure provides longitudinal confinement of the field in the nanofiber
guided modes. Due to the strong transverse confinement of the nanofiber
guided modes, the coupling between the atom and the guided modes will be
significantly enhanced and the strong-coupling regime can be realized, even
for moderate finesse \cite{cavity coupling}. Such PhC nanofiber system may
become a promising workbench for quantum non-linear optics and may serve as
nodes in quantum networks \cite{kimble}. Also the PhC nanofibers
may open up exciting new applications in lasing, optical switching and chemical/biological sensing.

In this paper, we present the fabrication of PhC nanofibers using
a femtosecond laser ablation technique. We demonstrate that thousands of periodic
nano-craters are fabricated on optical nanofibers by irradiating with just a
single femtosecond laser pulse. A key aspect of the fabrication is the
lensing effect of the nanofiber itself. We also demonstrate that such
periodic nano-craters on the nanofiber, induce strong modulation of refractive index
and act as a 1-D PhC.

\section{Experiments}

\subsection{Fabrication setup}

\begin{figure}[tbph]
\begin{center}
\includegraphics[height=8.5cm, width=8.5cm]{fig1.jpg}
\end{center}
\caption{(Color online) (a) The schematic diagram of the fabrication setup. The phase
mask splits the femtosecond laser beam into $\pm $ first orders which are
then recombined by the folding mirrors to create an interference pattern at
the nanofiber. A
cylindrical lens is used to line focus the femtosecond laser along the
nanofiber. A zero order block is used to avoid any residual zero order light
in the interference region. (b) The schematic diagram of the optical measurements.
The transmission and reflection spectra of the fabricated nanofiber samples are simultaneously measured by varying the
polarization of the input light (see text for details). NPBS
and OMA denote non-polarizing beam splitter and optical multi-channel analyzer, respectively.}
\label{fig1}
\end{figure}

A schematic diagram of the fabrication setup is shown in Fig. 1(a). The
femtosecond laser (Coherent Libra-HE) used for the fabrication has a central
wavelength ($\lambda $) of 400 nm and generates 120 fs pulses at a
repetition rate of 1 kHz with maximum pulse energy of 1.3 mJ. A Talbot
interferometer \cite{DUV}, consisting of a phase mask as beam splitter and
two folding mirrors, is built to create the two-beam interference pattern on
the nanofiber. The phase mask used for the fabrication has a uniform pitch ($%
\Lambda _{P}$) of 700 nm and was designed for zero order supression at 400
nm wavelength. The position of the folding mirrors is carefully chosen to
symmetrically recombine the first orders at the nanofiber position, thereby
creating an interference pattern with a period of $\Lambda _{P}/2=350$ nm.
A cylindrical lens is used to line-focus the femtosecond
laser beam along the nanofiber. The typical 1/e$^{2}$ beam size at the
nanofiber is 5.4 mm $\times $ 60 $\mu $m. The polarization of the
femtosecond laser beam is perpendicular to the nanofiber axis. Tapered
nanofibers with waist diameters ($2a$) of 450-650 nm are used for the
fabrication. After
the fabrication, the nanofiber samples are observed using a scanning electron
microscope (SEM).

\subsection{Alignment of the fabrication setup}
The interferometer alignment is highly sensitive to the optical path length
difference bewteen the two first orders. The path length difference must be
smaller than 36 $\mu $m as the spatial overlap is limited by the femtosecond
laser pulse width \cite{DUV, femtoreview1}. Hence, the positions of the
folding mirrors are controlled by micrometer stages and the mirror tilt
angles are precisely controlled. The alignment of the interferometer is
performed by observing the periodic ablation pattern on a glass plate, using
the SEM. After alignment of the interferometer, the
glass plate is replaced by the tapered nanofiber. For the fabrication, the
tapered fiber is mounted on a metallic holder having a central hole of 6 cm
length, so that the nanofiber is protected from contamination due to
ablation of the substrate. The overlap of the femtosecond laser beam with
the nanofiber is monitored by reducing the pulse energy to minimum value and
observing the scattered light from the nanofiber using a CCD camera. Maximum
overlap is achieved using a X-stage and a rotation-stage.

\subsection{Measurement of optical properties}
The optical properties of the
fabricated samples are characterized by measuring the transmission and
reflection spectra. A schematic of the
experimental setup is shown in Fig. 1(b). A supercontinuum light source (SuperK Extreme, NKT
Photonics) is launched into the tapered fiber sample and the spectrum of the
transmitted light is measured using a Fourier transform spectrometer
(FT-spectrometer). The resolution of the FT-spectrometer (Nicolet 8700,
Thermo Fisher Scientific) is $\sim 0.125$ cm$^{-1}$($\sim 0.01$ nm).
An inline fiber polarizer is used to control the polarization of the input
light. For confirmation, a part of the transmitted light is used for
polarization measurements. Simultaneously, the spectrum of the reflected
light is measured using an optical multi-channel analyzer (OMA). The
resolution of the OMA (QE65000, Ocean Optics) is $\sim 2$ nm. The
transmittance and reflectance values are calibrated using a CW Ti-Sapphire
laser source (MBR-110, Coherent Inc.).

\section{Results and Discussions}

\subsection{Nano-crater formation on nanofibers}

For the fabrication, we vary the pulse energy, repitition rate and
irradiation time (number of shots) to find the optimum conditions. It was found that for 1 kHz repitition rate, the
nanofiber can be completely destroyed within 1 second, even at pulse energy
of 200 $\mu $J. Even a repition rate of 100 Hz can still seriously damage the nanofiber.
However, we have found that by reducing the number shots to less than 20-shots, one could realize
controlled nanofabrication on the nanofiber. By systematically changing the number
of shots and pulse energy, we have found that single-shot irradiation is the
best condition to realize clean nanofabrication. In the following, we discuss the fabrication results
for multiple-shot and single-shot irradiation.

\subsubsection{Multiple-shot fabrication}

Figure 2(a) shows the SEM image of a
typical sample fabricated by 20-shot irradiation with a pulse energy of 370 $%
\mu $J. We have found that many ablated structures are formed on the nanofiber.
We must mention that these ablated structures are
formed, not on the irradiation side, but on the shadow side of the
nanofiber. As one can see the ablation pattern is quite irregular. However, a crucial observation is
that although the beam size along X-axis (60 $\mu $m) is much larger than
the nanofiber diameter, the ablation pattern is formed in a line along the
fiber axis (Z-axis).

Figure 2(b) shows the SEM image of a
typical sample fabricated by 3-shot irradiation with a pulse energy of 560 $%
\mu $J. One can see rather regular and periodic ablation pattern is formed on
the shadow side of the nanofiber.
However, the ablated structures are elliptical and are elongated along the Z-axis.
The typical size of an ablated structure is 90 nm $\times $ 210 nm. We must mention that
such ellipticity is observed throughout the ablation pattern, even for the weakest ablated structure.

\begin{figure}[tbph]
\begin{center}
\includegraphics[width=12.5cm]{fig2.jpg}
\end{center}
\caption{(Color online) SEM images of samples fabricated by multiple-shot irradiation.
(a) SEM image of a typical
sample fabricated by 20-shot irradiation. (b) SEM image of a typical
sample fabricated by 3-shot irradiation.
The ablated structures are observed on the shadow side of the nanofiber.}
\label{fig2}
\end{figure}

\subsubsection{Single-shot fabrication}

\begin{figure}[tbph]
\begin{center}
\includegraphics[width=12.5cm]{fig3.jpg}
\end{center}
\caption{(Color online) (a) SEM image of a typical sample
fabricated using single-shot irradiation. The inset shows the enlarged view.
The periodic nano-crater
structures are observed on the shadow side of the nanofiber. (b)
The cross-sectional image of a typical nano-crater measured by tilting the
nanofiber at an angle of 33$^{\circ }$.}
\label{fig2}
\end{figure}

Figure 3(a) shows the SEM image of a typical sample fabricated using
single-shot irradiation with a pulse energy of 630 $\mu $J. As one can
clearly see, periodic nano-crater structures are formed on the shadow side of the nanofiber.
Such periodic structures are quite systematically formed over a length of
few mm, consisting of thousands of such nano-craters. The inset shows the enlarged
view of the sample. The shape of the
nano-craters is almost circular and the diameter of a typical nano-crater is
around 210 nm. The nano-craters are formed with a periodicity of 350 nm,
which corresponds well to the interference fringe spacing. The depth of a
typical nano-crater is measured by cutting the nanofiber at the nano-crater
position and observing the cross-section using the SEM. The cross-sectional
view is shown in Fig. 3(b). The nano-crater has a bowl-like shape and the
depth is $\sim 120$ nm.

\subsubsection{Lensing effect of the nanofiber}

It should be noted that the ablated structures shown in Fig. 2 and Fig. 3 are
formed on the shadow side of the
nanofiber. We could not observe any particular structure on the irradiation
side. This suggests that the nanofiber itself acts as a cylindrical
lens and focuses the femtosecond laser beam on its shadow surface \cite%
{book1,trap}. Also as shown in Figs. 2(a) and (b), the ablation pattern is formed
exactly in a line along the fiber axis (Z-axis), confirming the lensing
effect of the nanofiber. The uniqueness of the present method is that the
lensing effect of the nanofiber makes it robust to any mechanical
instabilities in the transverse direction (X-axis). But
the instabilities along Z-axis can seriously affect the fabrication for multiple-shot irradiation, as the intensity pattern on the
nanofiber may differ for each shot. This is evident from the observed irregular ablation pattern for
20-shot irradiation (Fig. 2(a)), where the periodicity is completely washed out due to the instabilities along Z-axis.
Even for 3-shot irradiation (Fig. 2(b)) one can see that the ablated structures are elongated along the
fiber axis. However such instability does not affect the
fabrication in single-shot irradiation as the irradiation time is only 120
fs (i.e. pulse length). As a result periodic nanostructures with well
defined shape (shown in Fig. 3(a)) and periodicity are fabricated, without
taking any special care to suppress mechanical vibrations. Moreover, such single-shot irradiation
technique makes the fabrication method highly reproducible.
The lensing effect of the nanofiber may
explain the shape of the nano-craters. Although the field distribution along
the Z-axis is limited by the interference pattern, the circular shape of the
nano-craters is mainly due to focusing of the femtosecond laser along the
X-axis.

\subsection{Diameter distribution of nano-craters fabricated using single-shot irradiation}

\begin{figure}[tbph]
\begin{center}
\includegraphics[width=7.5cm]{fig4.jpg}
\end{center}
\caption{(Color online) The diameter distribution of the nano-craters for two nanofiber samples, (a) Sample 1 and (b) Sample 2,
fabricated using single-shot irradiation. The red squares denote the diameter of nano-craters
and the gray circles denote the corresponding diameter of the nanofiber.}
\label{fig3}
\end{figure}

The diameter distribution of the nano-craters for two nanofiber samples,
fabricated using single-shot irradiation, are
plotted in Figs. 4(a) and (b), along with the corresponding nanofiber diameter.
For Sample 1 the nano-craters are fabricated in the
tapered region, and the diameter of the nanofiber varies from 500 nm to 680
nm. A periodic array of nano-craters is formed in a region exceeding 2 mm
along the nanofiber axis. The diameter profile of the nano-craters shows a
peak-like behavior and the 1/e$^{2}$ width is $\sim 2.6$ mm. The diameter of
the nano-craters varies from 95 nm to 290 nm. On the other hand, for Sample
2 the fabrication is done near to the waist region, and the diameter of the
nanofiber is almost constant around 490 nm. The nano-craters are formed in a
region of 0.8 mm along the nanofiber. The diameter of the nano-craters
varies from 95 nm to 220 nm and the 1/e$^{2}$ width of the profile is $\sim
1.03$ mm. The observations also suggest rather smaller depth of nano-craters
for Sample 2 as compared to the Sample 1.

The observed 1/e$^{2}$ width of the diameter profile of the nano-craters is
smaller than the original beam profile. This is due to the fact that the
femtosecond laser induced ablation involves a multi-photon ionization
process \cite{femtoreview1, type2}. Hence the ablation process scales as a
higher power of the intensity and has some threshold condition. In order to understand the diameter
profile of the nano-craters one must also consider the nanofiber diameter, as the lensing effect of the
nanofiber depends on the size parameter ($2\pi a/\lambda $). In case of
Sample 2, the reduction in the length of the nano-crater array and the
smaller depths of the nano-craters may be understood from the decrease in
the lensing effect of the nanofiber due to thinner diameter \cite{book1,trap}. The observation of nano-crater diameters down to 95 nm suggests that
one can design various nanofabrications. The diameter distribution of the nano-craters
can be designed by controlling the intensity distribution and the size
parameter.

\subsection{Optical properties of nanofibers with periodic nano-craters}

\begin{figure}[tbph]
\begin{center}
\includegraphics[width=8.0cm]{fig5.jpg}
\end{center}
\caption{(Color online) Optical characteristics of nanofiber samples fabricated using
single-shot irradiation. (a) The transmission
spectrum of Sample 1. (b) The transmission and reflection spectra of Sample
2, measured for two orthogonal polarizations, X-polarization (black curves)
and Y-polarization (red curves).}
\label{fig4}
\end{figure}

Optical properties of the nanofiber samples, fabricated using
single-shot irradiation, are characterized by measuring the transmission and reflection spectra.
The transmission spectrum of Sample 1, shown in Fig. 5(a), shows a strong,
broad reflection band centered around 825 nm. The FWHM of the reflection
band is around 37 nm. In the wavelength region from 815 nm to 836 nm there
is almost no transmission, indicating the stop band. The reflectivity values
measured using a CW Ti-Sapphire laser confirms high reflectivity of $99\%$
in this stop band region. The transmission and reflection spectra of Sample
2 are shown in Fig. 5(b). The transmission and reflection spectra are
measured for two orthogonal linear polarizations, X-polarization (along
X-axis) and Y-polarization (along Y-axis). The transmission spectrum for the
X-polarization shows a reflection band centered around 823 nm and the FWHM
is 10 nm. The wavelength region from 819.5 nm to 824 nm denotes the stop
band. Note that the transmission spectrum for the Y-polarization is blue
shifted having a reflection peak at 817.5 nm and FWHM of 6.8 nm. The
spectrum is asymmetric and shows several small sharp peaks in the red side.
The reflection spectra almost match to the reflection band in the
transmission spectra except that the fine features are washed out due to the
resolution limit for the reflection measurements. The peak reflectivity
values are measured using a CW Ti-Sapphire laser. For the X-polarization, we
could observe almost 99\% reflection in the stop band region, whereas for the
Y-polarization the peak reflectivity is measured to be 80\%.

The observation of high reflectivity values in the stop band region, clearly demonstrates that the periodic nano-craters on the
nanofiber induce strong modulation of refractive index and act as a 1-D PhC.
The observed spectral width for Sample 1 is $3.7$ times broader than the
Sample 2. This is mainly due to the nanofiber diameter variation for Sample
1, as the diameter determines the propagation constant and the Bragg
resonance \cite{FIB 1, FIB 2}. In case of Sample 2, the observed blue shift
for the Y-polarization may be explained using the same argument \cite{FIB 1,
FIB 2}. Due to formation of nano-craters, the effective diameter along the Y-axis is reduced which results in such blue shift.
Also the transmission spectrum for the Y-polarization is asymmetric
and shows several small sharp peaks in the red side. Further investigations
are required to understand such fine spectral features.

The peak reflectivity value for the Y-polarization is measured to be $80\%$ whereas
the transmittance value at this wavelength is $12.5\%$. This suggests that,
for the Y-polarization, there is few percent scattering loss due to the
nano-craters \cite{FIB 1}. However, for the X-polarization, the scattering loss is negligible as is evident
from the high transmittance values away from the stop-band region. One should expect high
scattering loss as thousands of nano-craters are formed on the nanofiber.
The observed low scattering loss and high reflectivity is mainly due to the shape and shallow
depths of the nano-craters (shown in Fig. 3(b)). The observation of such excellent optical
properties, clearly suggests that one can design various types of 1-D PhC cavities on nanofibers using the present fabrication method.

\section{Conclusion}

In conclusion we have demonstrated the formation of periodic nano-craters on
a sub-wavelength diameter silica fiber using a femtosecond laser ablation
technique. The ablation is achieved by irradiating only a single femtosecond
pulse. Thanks to the lensing effect of the nanofiber, circular nano-crater
structures are formed on the shadow surface of the nanofiber. Such a
fabrication method may open new prospects in nanophotonics and
nanofabrication technologies. We have demonstrated that the periodic
nano-craters on the nanofiber act as a 1-D PhC. The PhC nanofiber\ system
can become a promising workbench for quantum non-linear optics and will open
new avenues in quantum information technology, by combining with
laser-cooled atoms or solid-state quantum emitters. Also the PhC nanofibers
may open up exciting new applications in lasing, optical switching and chemical/biological sensing.

\section*{Acknowledgements}
We are thankful to Mark Sadgrove and Makoto Morinaga for fruitful
discussions. Also we wish to thank Hitachi High-Tech Corp., Japan, for
helping in SEM measurements with SU8040. This work was supported by the
Japan Science and Technology Agency (JST) as one of the Strategic Innovation
projects.

\end{document}


\section{Introduction}

Let $(x_1, x_2, \ldots)$ be a positive sequence with finite sum $t = \sum_{i=1}^\infty x_i$. Its size-biased permutation (s.b.p) is the same sequence presented in a random order $(x_{\sigma_1}, x_{\sigma_2}, \ldots)$, where $\mathbb{P}(\sigma_1 = i) = \displaystyle \frac{x_i}{t}$, and for $k$ distinct indices $i_1, \ldots, i_k$,
\begin{equation}\mathbb{P}(\sigma_k = i_k| \sigma_1 = i_1, \ldots, \sigma_{k-1} = i_{k-1}) = \frac{x_{i_k}}{t - (x_{i_1} + \ldots + x_{i_{k-1}})}. \label{eqn:sbp} \end{equation}
An index $i$ with bigger `size' $x_i$ tends to appear earlier in the permutation, hence the name size-biased.
Size-biased permutation of a random sequence is defined by conditioning on the sequence values.

\subsection*{A brief history}
One of the earliest occurences of size-biased permutation is in social choice theory. Luce \cite{luce} studied distributions on the permutation on $n$ letters given by (\ref{eqn:sbp}) as a function of the $x_i$'s. These are the relative scores or desirabilities of the candidates, to be inferred through observing multiple rankings. This ranking model is now known as the Plackett-Luce model and it has wide applications \cite{plackett, plmeier, volkmer}.

Around the same time, biologists in population genetics were interested in inferring the distribution of alleles in a population through sampling. In these applications, $x_i$ is the abundance and $x_i/t$ is the relative abundance of the $i$-th species \cite{en78}. Size-biased permutation models the outcome of successive sampling, where one samples without replacement from the population and records the abundance of newly discovered species in the order that they appear. To account for the occurrence of new types of alleles through mutation and migration, they considered random abundance sequences and did not assume an upper limit to the number of possible types. Species sampling from random infinite sequence is sometimes known as \emph{size-biased random permutation}, a term coined by Patil and Taillie \cite{pt77}. The earliest work along this vein is perhaps that of McCloskey \cite{mc65}, who obtained results on the size-biased permutation of ranked jumps in a certain Poisson point process (p.p.p). The distribution of this ranked sequence is now known as the Poisson-Dirichlet $PD(0, \theta)$, and the distribution of its size-biased permutation is the $GEM(0, \theta)$; see Section \ref{subsec:parallel} for their definitions. This work was later generalized by Perman, Pitman and Yor \cite{pitmanyor}, who studied size-biased permutation of ranked jumps of a subordinator; see Section \ref{sec:sbp}.

The finite combinatorial version of size-biased sampling from $PD(0, \theta)$ is the Ewens sampling formula \cite{ew72}. Kingman \cite{ki78a} wrote: `One of the most striking results of recent theoretical research in population genetics is the sampling formula enunciated by Ewens and shown by (others) to hold for a number of different population models'. In studying this formula, Kingman initiated the theory of partition structure \cite{ki78a, ki78b}. Kingman showed that the Ewens sampling formula defines a particular partition structure by deletion of type; see \cite{ghp} for recent developments. Subsequent authors have studied partition structures and their representations in terms of exchangeable random partitions, random discrete distributions, random trees and associated random processes of fragmentation and coalescence, Bayesian statistics and machine learning. See \cite{csp} and references therein.
\subsection*{Organization}
This paper focuses on \emph{finite i.i.d size-biased permutation}, that is, the size-biased permutation of $n$ independent and identically distributed (i.i.d) random variables from some distribution $F$ on $(0, \infty)$. Our setting is a finite dimensional analogue of the size-biased permutation of ranked jumps of a subordinator studied in Perman-Pitman-Yor (PPY) \cite{pitmanyor}, as well as a special form of \emph{induced order statistics} \cite{david73, bhat74}; see Section \ref{sec:history} for a brief historical account of this field. We utilize this connection in our paper to arrive at new results. In Section \ref{sec:sbp} we study joint and marginal distribution of finite i.i.d size-biased permutation through a Markov chain. We draw connections between our settings and that of PPY \cite{pitmanyor}, and prove a converse of the stick-breaking representation of Patil and Taillie when $F$ is gamma (Proposition \ref{prop:pt}). Section \ref{sec:ios} we show that finite i.i.d size-biased permutation is a form of induced ordered statistics, and use this fact to derive previous distributional results. Comparisons to corresponding statements in Section \ref{sec:sbp} lead to a new beta-gamma identity when F is gamma (Corollary \ref{cor:newgb2}). As the sequence length tends to infinity, we derive asymptotics of the last $u$ fraction of finite i.i.d size-biased permutation in Section \ref{sec:asymp}, and that of the first few terms in Section \ref{sec:ppp}.

\subsection*{Notations} We shall write $gamma(a,\lambda)$ for a Gamma distribution whose density at $x$ is $\lambda^a x^{a-1}e^{-\lambda x}/\Gamma(a)$ for $x > 0$, and $beta(a,b)$ for the Beta distribution whose density at $x$ is $\displaystyle \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}x^{a-1}(1-x)^{b-1}$ for $x \in (0,1)$. For an ordered sequence $Y_n(k), k = 1, \ldots, n$, let $Y^{rev}_n(k) = Y_n(n-k+1)$ be the same sequence presented in reverse. For order statistics, we write $(Y^{\uparrow}(1), Y^{\uparrow}(2), \ldots)$ for the increasing sequence, and $(Y^{\downarrow}(1), Y^{\downarrow}(2), \ldots)$ for the decreasing sequence.  Throughout this paper we use $(X_1, \ldots, X_n)$ for the underlying i.i.d sequence, and $(X_n[1], \ldots, X_n[n])$ for its size-biased permutation.

\section{Connections to the work of Perman-Pitman-Yor\cite{pitmanyor}.}\label{sec:sbp}
PPY considered the size-biased permutation of the sequence of ranked jumps of a \emph{subordinator} $\tilde{T}$, that is, a non-decreasing process with no drift component, right continuous paths, stationary independent increments, $\tilde{T}(0) = 0$, and
$$ \tilde{T}(s) = \sum_i X_i \mathbf{1}_{\{\sigma_i \leq s\}} $$
for $ 0 \leq s \leq 1$, where $\{(\sigma_i, X_i), i \geq 1\}$ is the countable random set of points of a Poisson point process (p.p.p) on $(0, \infty)^2$ with intensity measure $ds \Lambda (dx)$. The $X_i$'s are the jumps of the subordinator $\tilde{T}$.
They assumed $\Lambda(dx) = \rho(x)\, dx$ for some density $\rho$, $\Lambda(0, \infty) = \infty, \Lambda(1, \infty) < \infty$, $\int_0^1 x \Lambda(dx) < \infty$. So $T > 0$ is an infinitely divisible random variable with L\'{e}vy measure $\Lambda$ and no drift component in its L\'{e}vy-Khitchine representation (see, for example, \cite[\S 15]{kallenberg}).

Proposition \ref{prop:converge} states that the PPY setup is the limit in distribution of finite i.i.d s.b.p. This is an application of the principle that convergence of order statistics is equivalent to convergence of size-biased permutation; see \cite{gnedin} for a precise statement.
\begin{prop}[Gnedin \cite{gnedin}]\label{prop:converge}
Let $X = (X^{\downarrow}(1), X^{\downarrow}(2), \ldots)$ be the infinite sequence of ranked jumps of a subordinator, whose sum $T = \sum_{i\geq 1}X^\downarrow(i)$ is a positive, infinitely divisible, a.s. finite random variable whose L\'evy-Khitchine representation has no drift component. Let $(X^n, n \geq 1)$ be an i.i.d positive triangular array, that is, $X^n = (X^n_1, \ldots, X^n_n)$ where $X^n_i$ are i.i.d, a.s. positive, and $T_n = \sum_{i=1}^nX^n_i \stackrel{d}{\rightarrow} T$ as $n \to \infty$. Then as $n \to \infty$,
$$ \mbox{s.b.p of } X^n \,\, \stackrel{f.d.d}{\longrightarrow} \mbox{ s.b.p of } X $$
meaning convergence in distributoin of finite dimensional distributions.
\end{prop}
\begin{proof} Recall that the sequence of decreasing order statistics $(X^{n, \downarrow}(1),\ldots, X^{n, \downarrow}(n))$ converges in distribution to $X$ \cite[\S 15]{kallenberg}. Since $T_n, T > 0$ a.s. and $T_n \stackrel{d}{\rightarrow} T$, the normalized sequence $(X^{n, \downarrow}(1)/T, \ldots, X^{n, \downarrow}(n)/T)$ converges in distribution to $X/T$. The result follows from
\cite[Theorem 3]{gnedin}.
\end{proof}

One can obtain another finite version of the PPY setup by setting $\Lambda(0, \infty) <~ \infty$, but this can be reduced to finite i.i.d s.b.p by conditioning. Specifically, $\tilde{T}$ is now a compound Poisson process, where the subordinator waits for an exponential time with rate $\Lambda(0, \infty)$ before making a jump $X$, whose length is independent of the waiting time and distributed as $\mathbb{P}(X \leq t) = \Lambda(0, t]/\Lambda(0, \infty)$ \cite{bertoin99}. If $(X_1, X_2, \ldots)$ is the sequence of successive jumps of $(\tilde{T}_s, s \geq 0)$, then $(X_1, X_2, \ldots, X_N)$ is the sequence of successive jumps of $(\tilde{T}_s, 0 \leq s \leq 1)$, where $N$ is a Poisson random variable with mean $\Lambda(0, \infty)$, independent of the jump sequence $(X_1, X_2, \ldots)$. For $N > 0$, properties of the size-biased permutation of $(X_1, \ldots, X_N)$ can be deduced from those of a finite i.i.d size-biased permutation by conditioning on $N$.

\subsection{Joint distribution, Markov property and stick-breaking}\label{subsec:joint.ppy}
Recall that $(X_1, \ldots, X_n)$ is an i.i.d sequence from distribution $F$, and $(X_n[1], \ldots, X_n[n])$ is its size-biased permutation. Let $T_{n-k} = X_n[k+1] + \ldots + X_n[n]$ denote the sum of the last $n-k$ terms in an i.i.d size-biased permutation of length $n$. We now derive the joint distribution of the first $k$ terms $X_n[1], \ldots, X_n[k]$ when $F$ has density $\nu_1$. The parallels in PPY are discussed in Section ~\ref{subsec:parallel}.

\begin{prop}[Barouch-Kaufman \cite{bkauf}] For $1 \leq k \leq n$, let $\nu_k$ be the density of $S_k$, the sum of $k$ i.i.d random variables with distribution $F$. Then
\begin{align}
& \mathbb{P}(X_n[1] \in dx_1, \ldots, X_n[k] \in dx_k) \nonumber \\
=& \frac{n!}{(n-k)!}\left(\prod_{j=1}^kx_j\, \nu_1(x_j)\, dx_j\right)\, \int_0^\infty \nu_{n-k}(s) \prod_{j=1}^k(x_j + \ldots + x_k+s)^{-1}\, ds. \label{eqn:joint.xnk}
\end{align}
\end{prop}
\begin{proof}
Let $\sigma$ denote the random permutation on $n$ letters defined by size-biased permutation as in (\ref{eqn:sbp}). Then there are $\frac{n!}{(n-k)!}$ distinct possible values for $(\sigma_1, \ldots, \sigma_k)$. By exchangability of the underlying i.i.d random variables $X_1, \ldots, X_n$, it is sufficient to consider $\sigma_1 = 1$, $\ldots$, $\sigma_k = ~k$. Note that
$$\mathbb{P}\left((X_1, \ldots, X_k) \in dx_1\, \ldots dx_k, \sum_{j=k+1}^nX_j \in ds\right) = \nu_{n-k}(s) \, ds \prod_{j=1}^k \nu_1(x_j) \, dx_j.$$
Thus, restricted to $\sigma_1 = 1$, $\ldots$, $\sigma_k = ~k$, the probability of observing $(X_n[1], \ldots, X_n[k]) \in dx_1\, \ldots dx_k$ and $T_{n-k} \in ds$ is precisely
$$\frac{x_1\, dx_1}{x_1 + \ldots + x_k + s}\, \frac{x_2\, dx_2}{x_2 + \ldots + x_k + s} \cdots \frac{x_k\, dx_k}{x_k + s}\, \nu_{n-k}(s)\, \left(\prod_{j=1}^k\nu_1(x_j)\, dx_j\right) ds. $$
By summing over $\frac{n!}{(n-k)!}$ possible values for $(\sigma_1, \ldots, \sigma_k)$, and integrating out the sum $T_{n-k}$, we arrive at (\ref{eqn:joint.xnk}).
\end{proof}
Note that $X_n[k] = T_{n-k+1} - T_{n-k}$ for $k = 1, \ldots, n-1$. Thus we can rewrite (\ref{eqn:joint.xnk}) in terms of the joint law of $(T_{n}, T_{n-1}, \ldots, T_{n-k})$:
\begin{equation}\label{eqn:joint.T}
\mathbb{P}(T_n \in dt_0, \ldots, T_{n-k} \in dt_k) = \frac{n!}{(n-k)!}\left(\prod_{i=0}^{k-1} \frac{t_i - t_{i+1}}{t_i}\nu_1(t_i - t_{i+1})\right)\, \nu_{n-k}(t_k)\, dt_0\ldots dt_k.
\end{equation}
Rearranging (\ref{eqn:joint.T}) yields the following result, which appeared as an exercise in \cite[\S 2.3]{cyor}.
\begin{cor}[Chaumont-Yor \cite{cyor}]
The sequence $(T_n, T_{n-1},\ldots, T_1)$ is an inhomogeneous Markov chain with transition probability
\begin{equation}\label{eqn:T.tp}
 \mathbb{P}(T_{n-k} \in ds| T_{n-k+1} = t) = (n-k+1)\, \frac{t-s}{t}\cdot\nu_1(t-s)\, \frac{\nu_{n-k}(s)}{\nu_{n-k+1}(t)}\, ds,
\end{equation}
for $k = 1, \ldots, n-1$. Together with $T_n \stackrel{d}{=} S_n$, equation (\ref{eqn:T.tp}) specifies the joint law in (\ref{eqn:joint.T}), and vice versa.
\end{cor}

An equivalent way to state (\ref{eqn:T.tp}) is that for $k \geq 1$, conditioned on $T_{n-k+1} = t$, $X_n[k]$ is distributed as the first size-biased pick out of $n-k+1$ i.i.d random variables conditioned to have sum $S_{n-k+1} = t$. This provides a recursive way to generate a finite i.i.d s.b.p: first generate $T_n$ (which is distributed as $S_n$). Conditioned on the value of $T_n$, generate $T_{n-1}$ via (\ref{eqn:T.tp}), let $X_n[1]$ be the difference. Now conditioned on the value of $T_{n-1}$, generate $T_{n-2}$ via (\ref{eqn:T.tp}), let $X_n[2]$ be the difference, and so on. Let us explore this recursion from a different angle by considering the ratio $\displaystyle W_{n,k} := \frac{X_n[k]}{T_{n-k+1}}$ and its complement, $\overline{W}_{n,k} = 1 - W_{n,k} = \displaystyle \frac{T_{n-k}}{T_{n-k+1}}$. For $k \geq 2$, note that
\begin{equation}\label{eqn:stickbreak}
 \frac{X_n[k]}{T_n} = \frac{X_n[k]}{T_{n-k+1}}\, \frac{T_{n-k+1}}{T_{n-k+2}}\cdots \frac{T_{n-1}}{T_n} = W_{n,k}\prod_{i = 1}^{k-1}\overline{W}_{n,i}.
\end{equation}
The variables $\overline{W}_{n,i}$ can be interpreted as residual fractions in a \emph{stick-breaking} scheme: start with a stick of length 1. Choose a point on the stick according to distribution $W_{n,1}$, `break' the stick into two pieces, discard the piece of length $W_{n,1}$ and rescale the remaining half to have length 1. Repeating this procedure $k$ times, and (\ref{eqn:stickbreak}) is the fraction broken off at step $k$ relative to the original stick length.

Together with $T_n \stackrel{d}{=} S_n$, one could use (\ref{eqn:stickbreak}) to compute the marginal distribution for $X_n[k]$ in terms of the ratios $\overline{W}_{n,i}$. In general the $W_{n,i}$ are not necessarily independent, and their joint distributions need to be worked out from (\ref{eqn:T.tp}). However, when $F$ has gamma distribution, $T_n, W_{n,1}, \ldots, W_{n,k}$ are independent, and (\ref{eqn:stickbreak}) leads to the following result of Patil and Taillie \cite{pt77}.

\begin{prop}[Patil-Taillie \cite{pt77}]\label{prop:pt} If $F$ has distribution $gamma(a, \lambda)$ for some $a, \lambda > 0$, then $T_n$ and the $W_{n,1}, \ldots, W_{n,n-1}$ in (\ref{eqn:stickbreak}) are mutually independent. In this case,
\begin{align*}
X_n[1] &= \gamma_0 \beta_1 \\
X_n[2] &= \gamma_0 \bar{\beta}_1 \beta_2 \\
& \ldots \\
X_n[n-1] &= \gamma_0 \bar{\beta}_1 \bar\beta_2\ldots \bar\beta_{n-2}\beta_{n-1} \\
X_n[n] &= \gamma_0 \bar{\beta}_1 \bar\beta_2\ldots \bar\beta_{n-1}
\end{align*}
where $\gamma_0$ has distribution $gamma(an,\lambda)$, $\beta_k$ has distribution $beta(a+1, (n-k)a)$, $\bar\beta_k = 1 - \beta_k$ for $1 \leq k \leq n-1$, and the random variables $\gamma_0, \beta_1, \ldots, \beta_n$ are independent.
\end{prop}
\begin{proof} By assumption $S_k$ has distribution $gamma(ak, \lambda)$. One substitutes the density of $gamma(ak, \lambda)$ for $\nu_k$ in (\ref{eqn:T.tp}), and the result follows by direct computation.
\end{proof}

In the subordinator setting, McCloskey \cite{mc65} proved the analogue of Proposition \ref{prop:pt}, and Perman-Pitman-Yor \cite{pitmanyor} proved the converse, see Proposition \ref{prop:mccloskey} below. Using the same idea, we obtain the following converse to Proposition \ref{prop:pt}, which appears to be new.
\begin{cor}\label{cor:pt} If $T_n$ is independent of $X_n[1]/T_n$, then $F$ is $gamma(a, \lambda)$ for some $a, \lambda > 0$.
\end{cor}
\begin{proof} Lukacs \cite{lukacs} proved that if $X, Y$ are non-degenerate, positive independent random variables, then $X+Y$ is independent of $\frac{X}{X+Y}$ if and only if both $X$ and $Y$ have gamma distributions with the same scale parameter. Note that
$$
\mathbb{P}(X_n[1]/T_n \in du, T_n \in dt) = n\, u\, \mathbb{P}\left(\frac{X_1}{X_1 + (X_2+\ldots+X_n)} \in du, T_n \in dt\right)
$$
and
$$\mathbb{P}(X_n[1]/T_n \in du) = n\, u\, \mathbb{P}\left(\frac{X_1}{X_1 + (X_2+\ldots+X_n)} \in du\right).$$
 Since $X_n[1]/T_n$ and $T_n$ are independent, $\frac{X_1}{X_1 + (X_2+\ldots+X_n)}$ is independent of $T_n$. The conclusion follows by applying Lukacs' theorem to the pair $X_1$ and $(X_2 + \ldots + X_n)$.
\end{proof}

\subsection{Parallels in PPY \cite{pitmanyor}} \label{subsec:parallel}
In the subordinator setting of PPY, let $\tilde{T}_k$ denote the remaining sum after removing the first $k$ terms of the size-biased permutation of the infinite sequence $(X^\downarrow(1), X^\downarrow(2), \ldots)$. The sequence $(\tilde{T}_0, \tilde{T}_1, \ldots)$ is a Markov chain with stationary transition probabilities \cite[equation 2.a]{pitmanyor}
$$
 \mathbb{P}(\tilde{T}_{1} \in dt_1| \tilde{T}_0 = t) = \frac{t-t_1}{t}\cdot\rho(t-t_1)\, \frac{\nu(t_1)}{\nu(t)}\, dt_1.
 $$

Conditionally given $\tilde{T}_0 = t_0, \tilde{T}_1 = t_1, \ldots, \tilde{T}_n = t_n$, the sequence of remaining terms in the s.b.p $(X[n+1], X[n+2], \ldots)$ is distributed as $(X^\downarrow(1), X^\downarrow(2), \ldots, )$ conditioned on $\sum_{i\geq 1}X^\downarrow(i) = t_n$, independent of the first $n$ size-biased picks \cite[Theorem 4.2]{pitmanyor}. The stick-breaking representation in (\ref{eqn:stickbreak}) now takes the form
\begin{equation}\label{eqn:stick2}
 \frac{X[k]}{\tilde{T}_0} = W_{k}\prod_{i=1}^{k-1}\overline{W}_{i},
\end{equation}
where $X[k]$ is the $k$-th size-biased pick, and $W_{i} = \frac{X[i]}{\tilde{T}_{i-1}}$, $\overline{W}_{i} = 1 - W_{i} = \frac{\tilde{T}_i}{\tilde{T}_{i-1}}$. Proposition ~\ref{prop:pt} and Corollary \ref{cor:pt} parallel the following result.
\begin{prop}[McCloskey \cite{mc65} and PPY \cite{pitmanyor}]\label{prop:mccloskey}
The random variables $\tilde{T}_0$ and $W_{1}, W_2, \ldots$ in (\ref{eqn:stick2}) are mutually independent if and only if $\tilde{T}_0$ has distribution $gamma(a, \lambda)$ for some $a, \lambda > 0$. In this case, the $W_{i}$ are i.i.d with distribution $beta(1, a)$ for $i = 1, 2, \ldots$.
\end{prop}

We take a small detour to explain some results related to Proposition \ref{prop:pt} and \ref{prop:mccloskey}
on characterization of size-biased permutations. For a random discrete distribution prescribed by its probability mass function $(P_k) = (P_1, P_2, \ldots)$ with $\sum_iP_i = 1$, $P_i \geq 0$, let $(\tilde{P}_k)$ be the s.b.p. of $(P_k)$. One may ask when is a given distribution $(Q_k)$ the s.b.p of some distribution $(P_k)$. Clearly $(\widetilde{\tilde{P}}_k) \stackrel{d}{=} (\tilde{P}_k)$ for any $(P_k)$, thus this question is equivalent to characterizing random discrete distributions on $\mathbb{N}$ which are invariant under size-biased permutation (ISBP).
Pitman \cite[Theorem 4]{pitman96} gave a complete answer in terms of symmetry of a certain function of the finite dimensional distributions. Furthermore, Pitman proved a complete characterization of ISBP when $P_k$ can be written as the right hand side of (\ref{eqn:stick2}) with $W_{1}, W_2, \ldots$ independent. In this situation, apart from some limiting cases, $(P_k)$ is ISBP if and only if $W_{i}$ is distributed as $beta(1-\alpha, \theta + i\alpha)$, $i = 1, 2, \ldots$ for certain pairs of real numbers $(\alpha, \theta)$. The two main cases are $(0 \leq \alpha < 1, \theta > -\alpha)$ and $(\alpha = -a < 0, \theta = na)$ for some $n = 1, 2, \ldots$. In both settings, $(P_k)$ is known as the $GEM(\alpha, \theta)$ distribution. The abbreviation GEM was introduced by Ewens, which stands for Griffiths-Engen-McCloskey. The McCloskey case of Proposition \ref{prop:mccloskey} is $GEM(0, \theta)$, and the Patil-Taillie case of Proposition \ref{prop:pt} is $GEM(-a, na)$. The sequence obtained by ranking a $GEM(\alpha, \theta)$ sequence is called a Poisson-Dirichlet distribution with parameters $(\alpha, \theta)$ \cite{pitmanyor}. The $GEM$ distribution has a generative description known as the Chinese restaurant process \cite[\S 3]{csp} and has applications in Bayesian statistics and machine learning, see \cite{csp}.

\section{Connections to induced order statistics}\label{sec:ios}
When $n$ i.i.d pairs $(X_i, Y_i)$ are ordered by their $Y$-values, the corresponding $X_i$'s are called the \emph{induced order statistics} of the vector $Y$, or its \emph{concomitants}. Gordon \cite{gordon83} first proved the following result for finite $n$ which shows that finite i.i.d s.b.p is a form of induced order statistics. Here we state the infinite sequence version, which is a special case of PPY\cite[Lemma 4.4]{pitmanyor}.

\begin{prop}[PPY \cite{pitmanyor}]\label{prop:coupling} Let $(x_1, x_2, \ldots)$ be a fixed positive sequence with finite sum $t = \sum_{i=1}^\infty x_i$. For $i = 1, 2,\ldots$, let $Y_i = \epsilon_i/x_i$ where $\epsilon_i$ are i.i.d standard exponentials. Let $(Y^{\uparrow}(1),Y^{\uparrow}(2), \ldots)$ be the increasing order statistics of the $Y_i$'s, and let $X^\ast(k)$ be the value of the $x_i$ such that $Y_i$ is $Y^{\uparrow}(k)$. Then $(X^\ast(k), k = 1, 2, \ldots)$ is a size-biased permutation of $(x_1, x_2, \ldots)$. In particular, the size-biased permutation of a positive i.i.d sequence $(X_1, \ldots, X_n)$ is distributed as the induced order statistics of the sequence $(Y_i = \epsilon_i/{X_i}, 1 \leq i \leq n)$ for an independent sequence of i.i.d standard exponentials $(\epsilon_1, \ldots, \epsilon_n)$.
\end{prop}
\begin{proof} Note that the $Y_i$'s are independent exponentials with rates $x_i$. Let $\sigma$ be the random permutation such $Y_{\sigma(i)} = Y^\uparrow(i)$. Note that $X^\ast(k) = x_{\sigma(k)}$. Then
$$\mathbb{P}(\sigma(1) = i) = \mathbb{P}(Y_i = \min\{Y_j, j = 1, 2, \ldots\}) = \frac{x_i}{t},$$
thus $X_n^\ast(1) \stackrel{d}{=} x[1]$. In general, for distinct indices $i_1, \ldots, i_k$, by the memoryless property of the exponential distribution,
\begin{align*}
& \, \mathbb{P}(\sigma(k) = i_k  | \sigma(1) = i_1, \ldots, \sigma(k) = i_{k-1} ) \\
=& \, \mathbb{P}\left(Y_{i_k} = \min\{Y_{\sigma(j)}, j \geq k \} | \sigma(1) = i_1, \ldots, \sigma(k) = i_{k-1}\right) = \frac{x_{i_k}}{t - \sum_{j=1}^{k-1} x_{i_j}}.
\end{align*}
Induction on $k$ completes the proof.
\end{proof}

This proposition readily supplies simple proofs for joint, marginal and asymptotic distributions of i.i.d s.b.p, as explicitly computed in this section. For instance, the proof of the following nesting property of i.i.d s.b.p, which can be cumbersome, amounts to i.i.d thinning.
\begin{cor}Consider a finite i.i.d s.b.p $(X_n[1], \ldots, X_n[n])$ from a distribution $F$. For $1 \leq m \leq n$, select $m$ integers $a_1 < \ldots < a_m$ by uniform sampling from $\{1, \ldots, n\}$ without replacement. Then the subsequence $\{X_n[a_j], 1 \leq j \leq m\}$ is jointly distributed as a finite i.i.d s.b.p of length ~$m$ from $F$.
\end{cor}

\subsection{Joint and marginal distribution revisited} \label{subsec:ios.joint}

Proposition \ref{prop:coupling} immediately yields the joint distribution of the $X_n[k]$'s, which is a different formula to (\ref{eqn:joint.xnk}).
\begin{prop}\label{lem:joint.x.u}
$(X_n[k], k = 1, \ldots, n)$ is distributed as the first coordinate of the sequence of pairs $((X^{\ast}_n(k), U_n^\downarrow(k)), k =1, \ldots, n)$, where $U_n^\downarrow(1) \geq \ldots \geq U^\downarrow_n(n)$ is a sequence of uniform order statistics, and conditional on $(U^\downarrow_n(k) = u_k, 1 \leq k \leq n)$, the $X^{\ast}_n(k)$ are independent with distribution $(G_{u_k}(\cdot), k = 1, \ldots, n)$, where
\begin{equation} \label{eqn:Gu}
 G_u(dx) = \frac{xe^{-\phi^{-1}(u)x}\, F(dx)}{-\phi'(\phi^{-1}(u))}.
\end{equation}
Here $\phi$ is the Laplace transform of $X$, that is, $\phi(y) =~ \int_0^\infty e^{-yx}\, F(dx)$, $\phi'$ its derivative and $\phi^{-1}$ its inverse function.
\end{prop}
\begin{proof} Let $(X_1, \ldots, X_n)$ be $n$ i.i.d draws from $F$, $Y_i = \epsilon_i/X_i$ for $n$ i.i.d standard exponentials $\epsilon_i$, independent of the $X_i$'s. Note that $\{(X_i, Y_i), 1 \leq i \leq n\}$ is an i.i.d sample from the joint distribution $ F(dx)[xe^{-yx}\, dy].$ Thus $Y_i$ has marginal density
\begin{equation}
P(Y_i \in dy) = -\phi'(y)\, dy, \hspace{1em} 0 < y < \infty, \label{1.b}
\end{equation}
and its distribution function is $F_Y = 1 - \phi$.  Given $\{Y_i =  y_i, 1 \leq i \leq n\}$, the $X_n^\ast(i)$ defined in Proposition \ref{prop:coupling} are independent with conditional distribution $\widetilde{G}(y_i, \cdot)$ where
\begin{equation}
\widetilde{G}(y, dx) = \frac{xe^{-yx}\, F(dx)}{-\phi'(y)}. \label{eqn:Gy}
 \end{equation}

Equation (\ref{eqn:Gu}) follows from writing the order statistics as the inverse transforms of ordered uniform variables
\begin{equation}\label{eqn:ynun}
(Y_n^\uparrow(1), \ldots, Y^\uparrow_n(n)) \stackrel{d}{=} (F_Y^{-1}(U^\downarrow_n(n)), \ldots, F_Y^{-1}(U^\downarrow_n(1))) \stackrel{d}{=} (\phi^{-1}(U^\downarrow_n(1)), \ldots, \phi^{-1}(U^\downarrow_n(n)))
\end{equation}
where $(U^\downarrow_n(k), k = 1, \ldots, n)$ is an independent decreasing sequence of uniform order statistics. Note that the minus sign in (\ref{1.b}) results in the reversal of the sequence $U_n$ in the second equality of (\ref{eqn:ynun}).
\end{proof}
\begin{cor} \label{cor:marginal}
For $1 \leq k \leq n$ and $0 < u < 1$, let
\begin{equation}\label{2.a}
f_{n,k}(u) = \frac{\mathbb{P}(U_{n}(k) \in du)}{du} = n\binom{n-1}{k-1}u^{n-k}(1-u)^{k-1}
\end{equation}
be the density of the $k$-th largest of the $n$ uniform order statistics $(U_n(i), i = 1, \ldots, n)$. Then
\begin{equation}\label{2.c}
\frac{\mathbb{P}(X_n[k] \in dx)}{xF(dx)} = \int_0^\infty e^{-xy}f_{n,k}(\phi(y))\, dy.
\end{equation}
\end{cor}
In particular, for the first and last values,
\begin{align*}
\frac{\mathbb{P}(X_n[1] \in dx)}{xF(dx)} &= n\int_0^\infty e^{-xy}\phi(y)^{n-1}\, dy \\
\frac{\mathbb{P}(X_n[n] \in dx)}{xF(dx)} &= n\int_0^\infty e^{-xy}(1-\phi(y))^{n-1}\, dy.
\end{align*}

\begin{ex}\label{ex:gamma1}
Suppose $F$ is $gamma(a,1)$. Then $\phi(y) = (\frac{1}{1+y})^a$, and $\phi^{-1}(u) =~ u^{-1/a} - 1.$
Hence $G_u$ in (\ref{eqn:Gu}) is
$$G_u(dx) = \frac{x}{au^{(a+1)/a}}\, e^{-(u^{-1/a} - 1)x}\, F(dx) = \frac{x^a}{\Gamma(a+1)}u^{-(a+1)/a}e^{-xu^{-1/a}}.$$
That is, $G_u$ is $gamma(a+1, u^{-1/a})$.
\end{ex}

\subsection{A new beta-gamma identity}\label{subsec:ios.gamma}
When $F$ is $gamma(a,\lambda)$, Lemma \ref{lem:joint.x.u} gives the following result, which is a remarkable complement to the Patil-Taillie representation in Proposition \ref{prop:pt}.

\begin{prop}\label{prop:newgb}
Suppose $F$ is $gamma(a,\lambda)$. Then $G_u$ is $gamma(a+1, \lambda u^{-1/a})$, and
\begin{equation}\label{eqn:xrev.g}
 (X_n[k], k = 1, \ldots, n) \stackrel{d}{=} ([U^\downarrow_n(k)]^{1/a}\gamma_k, k = 1, \ldots, n)
 \end{equation}
where $\gamma_1, \ldots, \gamma_n$ are i.i.d $gamma(a+1, \lambda)$ random variables, independent of the sequence of decreasing uniform order statistics $(U^\downarrow_n(1), \ldots, U^\downarrow_n(n))$. Alternatively, jointly for $k =~ 1, \ldots, n$
\begin{align*}
X_n^{rev}[1] &= \gamma_1 \beta_{an, 1} \\
X_n^{rev}[2] &= \gamma_2 \beta_{an, 1}\beta_{an-a, 1} \\
& \ldots \\
X_n^{rev}[n-1] &= \gamma_{n-1} \beta_{an, 1}\beta_{an-a, 1}\cdots \beta_{2a, 1} \\
X_n^{rev}[n] &= \gamma_{n} \beta_{an, 1}\beta_{an-a, 1}\cdots \beta_{a, 1}.
\end{align*}
where the $\beta_{an - ia,1}$ for $i = 0, \ldots, n-1$ are distributed as $beta(an-ia, 1)$, and they are independent of each other and the $\gamma_k$'s.
\end{prop}
\begin{proof} The distribution $G_u$ is computed in the same way as in Example \ref{ex:gamma1} and (\ref{eqn:xrev.g}) follows readily from Proposition \ref{lem:joint.x.u}.
\end{proof}
A direct comparison of the two different representations in Proposition \ref{prop:pt} and \ref{prop:newgb} creates $n$ distributional identities. For example, the equality $X_n[1] = X^{rev}_n[n]$ shows that the following two means of creating a product of independent random variables produce the same result in law:
\begin{equation}\label{eqn:beta.ini}
 \beta_{a+1, (n-1)a} \, \gamma_{an, \lambda} \stackrel{d}{=} \beta_{an, 1} \, \gamma_{a+1, \lambda}
 \end{equation}
 where $\gamma_{r,\lambda}$ and $\beta_{a,b}$ denote random variables with distributions $gamma(r, \lambda)$ and $beta(a,b)$, respectively. In fact, by comparing the consecutive ratios $\frac{X_n[i+1]}{X_n[i]}$ for $i = 1, \ldots, n-1$, the other $n-1$ identities all reduce to the following single equation after setting $ia - a = b$.
\begin{cor}\label{cor:newgb2} For $a, b, \lambda > 0$,
\begin{equation}\label{eqn:beta.ratio}
\frac{1 - \beta_{a+1,a+b}}{\beta_{a+1,a+b}}\beta_{a+1, b} \stackrel{d}{=} \frac{\gamma_{a+1,\lambda}}{\gamma'_{a+1,\lambda}}\beta_{a+b,1},
 \end{equation}
where the random variables $\beta_{a+1,a+b}, \beta_{a+1,b}, \beta_{a+b,1}, \gamma_{a+1,\lambda}, \gamma'_{a+1,\lambda}$ are mutually independent, beta and gamma distributed with parameters indicated by subscripts.
\end{cor}
The validity of this identity can be checked by comparing moments. Conversely, (\ref{eqn:beta.ini}) and (\ref{eqn:beta.ratio}) allow one to go between Proposition \ref{prop:newgb} and the Proposition \ref{prop:pt}.

\section{Asymptotics of the last $u$ fraction of the size-biased permutation}\label{sec:asymp}
In this section we derive Glivenko-Cantelli and Donsker-type theorems for the distribution of the last $u$ fraction of terms in a finite i.i.d s.b.p. These are special cases of more general statements which hold for arbitrary induced order statistics in $d$ dimensions (see Section \ref{sec:history}).
Features pertaining to i.i.d s.b.p are presented in the following lemma, which has an interesting successive sampling interpretation; see Section \ref{subsec:heuristic}.
\begin{lem}\label{lem:ode} Suppose $F$ has support on $[0, \infty)$ and finite mean. For $u \in (0,1)$, define
\begin{equation}\label{eqn:Fu}
F_u(dx) = \frac{e^{-x\phi^{-1}(u)}}{u}\, F(dx)
\end{equation}
and extend the definition to $\{0, 1\}$ by continuity, where $\phi$ is the Laplace transform of $F$ as in Proposition \ref{lem:joint.x.u}. Then $F_u$ is a probability distribution on $[0, \infty)$ for all $u \in [0,1]$, and $G_u$ in (\ref{eqn:Gu}) satisfies
\begin{equation}\label{eqn:Gu.is.sb}
G_u(dx) = xF_u(dx)/ \mu_u
\end{equation}
where $\mu_u = \int x F_u(dx) = \frac{-\phi'(\phi^{-1}(u))}{u}$. Furthermore,
\begin{equation}\label{eqn:integral} \int_0^u G_s(dx)\, ds = F_u(dx) \end{equation}
for all $s \in [0,1]$. In other words, the density
$$f(u,x) = F_u(dx) / F(dx) = u^{-1}e^{-x\phi^{-1}(u)}$$
of $F_u$ with respect to $F$ solves the differential equation
\begin{equation}\label{eqn:evo}
 \frac{d}{du}[uf(u,x)] = \frac{-xf(u,x)}{\mu_u} \end{equation}
with boundary condition $f(1,x) \equiv 1$.
\end{lem}
\begin{proof} By direct computation.
\end{proof}

We now state a Glivenko-Cantelli-type theorem which applies to size-biased permutation of finite deterministic sequences. Versions of this result are known in the literature \cite{bnw, holst, gordon83, rosen}, see discussions in Section \ref{sec:history}. We offer an alternative proof using induced order statistics.

\begin{thm}\label{thm:gc} Let $((x_i^n, 1 \leq i \leq n), 1 \leq n)$ be a deterministic triangular array of positive numbers  with corresponding c.d.f sequence $(E_n, 1 \leq n)$. Suppose
\begin{equation}
\sup_x|E_n(x) - F(x)| \to 0 \mbox{ as } n \to \infty \label{eqn:sup.condition}
\end{equation}
for some distribution $F$ on $(0, \infty)$. Let $u \in (0, 1]$. Let $E_{n,u}(\cdot)$ be the empirical distribution of the first $\lfloor nu \rfloor$ terms in a size-biased permutation of the sequence $(x_i^n, 1 \leq i \leq n)$. Then for each $\delta \in (0,1)$,
\begin{equation}\label{eqn:gc.main}
\sup_{u \in [\delta, 1]} \sup_\mathbf{I} | E_{n,u}(\mathbf{I}) - F_u(\mathbf{I})| \stackrel{a.s.}{\longrightarrow} 0 \mbox{ as } n \to \infty,
\end{equation}
where $\mathbf{I}$ ranges over all subintervals of $(0, \infty)$.
\end{thm}

\begin{proof} Define $Y^n_i = \epsilon_i/x^n_i$ for $i = 1, \ldots, n$ where $\epsilon_i$'s are i.i.d standard exponentials as in Proposition \ref{prop:coupling}. Let $H_n$ be the empirical distribution function (e.d.f) of the $Y^n_i$'s,
$$H_n(y) := \frac{1}{n}\sum_{i=1}^n\mathbf{1}_{\{Y^n_i < y\}}.$$
Let $J_n$ denote the e.d.f of $(x^n_i, Y^n_i)$. By Proposition \ref{lem:joint.x.u},
\begin{equation}\label{eqn:eemp}
E_{n,u}(\mathbf{I}) = \frac{n}{\lfloor nu \rfloor}\, \frac{1}{n} \sum_{i=1}^n\mathbf{1}_{\{x^n_i \in \mathbf{I}\}}\mathbf{1}_{\{Y^n_i < H_n^{-1}(1-u)\}} = \frac{n}{\lfloor nu \rfloor}J_n(\mathbf{I} \times [0, H_n^{-1}(1-u)]).
\end{equation}
Fix $\delta \in (0, 1)$, and let $u \in [\delta, 1]$. Let $\phi$ be the Laplace transform of $F$ and $J$ the joint law of $(X, \epsilon/X)$, where $X$ is a random variable with distribution $F$, and $\epsilon$ is an independent standard exponential. Note that $\frac{1}{u}J(\mathbf{I} \times [0, \phi^{-1}(u)]) = F_u(\mathbf{I}).$
Thus
\begin{align}
 \mathbb{E}_{n, u}(\mathbf{I}) - F_u(\mathbf{I}) = & \left(\frac{n}{\lfloor nu \rfloor}J_n(\mathbf{I} \times [0, H_n^{-1}(1-u)]) - \frac{n}{\lfloor nu \rfloor}J_n(\mathbf{I} \times [0, \phi^{-1}(u)])\right) \nonumber \\
& + \left(\frac{n}{\lfloor nu \rfloor}J_n(\mathbf{I} \times [0, \phi^{-1}(u)]) - \frac{1}{u}J(\mathbf{I} \times [0, \phi^{-1}(u)])\right). \label{eqn:split}
\end{align}
Let us consider the second term. Note that
\begin{align*}
J_n(\mathbf{I} \times [0, \phi^{-1}(u)]) &= \int_0^\infty e^{-t\phi^{-1}(u)}\mathbf{1}_{\{t \in \mathbf{I}\}}\, E_n(dt).
\end{align*}
Since $E_n$ converges to $F$ uniformly and $e^{-t\phi^{-1}(u)}$ is bounded for all $t \in (0, \infty)$ and $u \in [\delta, 1]$,
$$\sup_{u \in [\delta, 1]}\sup_\mathbf{I}\left|\frac{n}{\lfloor nu \rfloor}J_n(\mathbf{I} \times [0, \phi^{-1}(u)]) - \frac{1}{u}J(\mathbf{I} \times [0, \phi^{-1}(u)])\right| \stackrel{a.s.}{\longrightarrow} 0 \mbox{ as } n \to \infty.$$
Let us consider the first term.
Since $J_n$ is continuous in the second variable, it is sufficient to show that
\begin{equation}\label{eqn:hncont}
 \sup_{u \in [\delta, 1]}|H_n^{-1}(1-u) - \phi^{-1}(u)| \stackrel{a.s.}{\longrightarrow} 0 \mbox{ as } n \to \infty.
 \end{equation}
 To achieve this, let $A_n$ denote the `average' measure
$$A_n(y) := \frac{1}{n}\sum_{i=1}^n\mathbb{P}(Y^n_i < y) = 1 - \int_0^\infty e^{-xy}\, dE_n(x).$$
A theorem of Wellner \cite[Theorem 1]{wellner} states that if the sequence of measures $(A_n, n \geq 1)$ is tight, then the Prohorov distance between $H_n$ and $A_n$ converges a.s. to $0$ and $n \to \infty$. In this case, since $E_n$ converges to $F$ uniformly, $A_n$ converges uniformly to $1 - \phi$. Thus $H_n$ converges uniformly to $1-\phi$, and (\ref{eqn:hncont}) follows.
\end{proof}

In particular, when $E_n$ is the e.d.f of $n$ i.i.d picks from $F$, then (\ref{eqn:sup.condition}) is satisfied a.s. by the Glivenko-Cantelli theorem. Thus Theorem \ref{thm:gc} implies that $F_u$ is the limit of the empirical distribution function (e.d.f) of the last $u$ fraction in a finite i.i.d s.b.p.
\begin{cor}\label{cor:emp}For $u \in (0,1]$, let $F_{n,u}(\cdot)$ denote the empirical distribution of the last $\lfloor nu \rfloor$ values of an i.i.d s.b.p. with length $n$. For each $\delta \in (0,1)$, as $n \to \infty$,
\begin{equation}\label{emp.u}
\sup_{u \in [\delta, 1]} \sup_\mathbf{I} | F_{n,u}(\mathbf{I}) - F_u(\mathbf{I})| \stackrel{a.s.}{\longrightarrow} 0,
\end{equation}
where $\mathbf{I}$ ranges over all subintervals of $(0, \infty)$.
\end{cor}

Under the settings in Theorem \ref{thm:gc}, one can also derive a Donsker-type result via entropy methods described in \cite[\S 2]{vaart2}, complementing similar results in the literature \cite[\S 5]{bnw}. In particular, the first term of (\ref{eqn:split}) depends on the flunctuation of empirical uniform order statistics around its limit. Its scaling depends on the rate of convergence of $E_n$ to $F$. The second term of (\ref{eqn:split}) converges in distribution after rescaled by $\sqrt{n}$, by a version of Donsker's theorem, to a Brownian bridge. We omit the details. A similar decomposition applies in Theorem ~\ref{thm:fclt}.

\subsection{A heuristic interpretation} \label{subsec:heuristic}
Since $X_n^{rev}[\lfloor u \, n \rfloor]$ converges in distribution to $G_u$ for $u \in [0, 1]$, Corollary \ref{cor:emp} lends a sampling interpretation to Lemma \ref{lem:ode}. Equation (\ref{eqn:evo}) has the heuristic interpretation as characterizing the evolution of the mass at $x$ over time $u$ in a successive sampling scheme. To be specific, consider a successive sampling scheme on a large population of $N$ individuals, with species size distribution $H$. Scale time such that at time $u$, for $0 \leq u \leq 1$, there are $Nu$ individuals (from various species) remaining to be sampled. Let $H_u$ denote the distribution of species sizes at time $u$, and fix the bin $(x, x + dx)$ of width $dx$ on $(0, \infty)$. Then $NuH_u(dx)$ is the number of individuals whose species size lie in the range $(x, x + dx)$ at time $u$. Thus $\frac{d}{du}NuH_u(dx)$ is the rate of individuals to be sampled from this range of species size at time $u$. The probability of an individual whose species size is in $(x, x+dx)$ being sampled at time $u$ is $\frac{xH_u(dx)}{\int_0^\infty xH_u(dx)}$. As we scaled time such that $u \in [0,1]$, in time $du$ we sample $N \, du$ individuals. Thus
$$ \frac{d}{du}NuH_u(dx) = -N\, \frac{xH_u(dx)}{\int_0^\infty xH_u(dx) }.$$
Let $f(u,x) = H_u(dx) / H(dx)$, then as a function in $u$, the above equation reduces to (\ref{eqn:evo}).

\subsection{Functional central limit theorem}\label{subsec:fclt}
We now state a functional central limit theorem for i.i.d s.b.p. This is a special case of the general result of Davydov and Egorov \cite{egorov89} for induced order statistics.

\begin{thm}[Davydov and Egorov, \cite{egorov89}]\label{thm:fclt} Suppose the first two moments of $F$ are finite. For a distribution $H$, let $\mu(H)$, $\sigma(H)$ denote its mean and standard deviation. For $u \in (0, 1]$, define
\begin{align*}
\xi_n(u) &= \frac{1}{n}\sum_{j=1}^{\lfloor nu \rfloor} X^{rev}_n[j] \\
\eta_n(u) &= \frac{1}{n}\sum_{j=1}^n X_j \mathbf{1}(Y_j < \phi^{-1}(u)) \\
m(u) &= \int_0^u \mu(G_s) \, ds = \mu(F_u).
\end{align*}
Then $\xi_n \to m$, $\eta_n \to m$ uniformly on $[0,1]$, and
$$\sqrt{n}(\xi_n - m) \Rightarrow \alpha, \hspace{2em} \sqrt{n}\, (\eta_n - m) \Rightarrow \eta $$
in the Skorokhod topology, where
$$\alpha(u) = \int_0^u \sigma(G_s) dW(s), \hspace{1em} \eta(u) = \alpha(u) + \int_0^u m(s) \, dV(s).$$
$W$ is a standard Brownian motion, and $V$ is a Brownian bridge, independent of $W$.
\end{thm}
The difference in $\xi_n$ and $\eta_n$ is the fluctuation of empirical uniform order statistics around its limit. Recall the analogue discussed at the end of Theorem \ref{thm:gc}. The proof of Theorem \ref{thm:fclt} can be found in \cite{egorov89}, together with similar statements on functional law of iterated logarithm for the processes $\eta_n, \xi_n$.
\subsection{Historical notes on induced order statistics and successive sampling} \label{sec:history}
Induced order statistics were first introduced by David \cite{david73} and independently by Bhattacharya \cite{bhat74}. Typical applications stem from modeling an indirect ranking procedure, where subjects are ranked based on their $Y$-attributes although the real interest lies in ranking their $X$-attributes, which are difficult to obtain at the moment where the ranking is required.\footnote{One often uses $X$ for the variable to be ordered, and $Y$ for the induced variable, with the idea that $Y$ is to be predicted. Here we use $X$ for the induced order statistics since $X_n[k]$ has been used for the size-biased permutation. The role of $X$ and $Y$ in our case is interchangeable, as evident when one writes $X_iY_i = \epsilon_i$.} For example in cattle selection, $Y$ may represent the genetic makeup, for which the cattle are selected for breeding, and $X$ represents the milk yields of their female offspring. Thus a portion of this literature focuses on comparing distribution of induced order statistics to that of usual order statistics \cite{dayang77, yang77, nagada94, hall94}.
The most general statement on asymptotic distributions is obtained by Davydov and Egorov \cite{egorov89}, who proved functional central limit theorem and functional law of the iterated logarithm for the process $S_{n,u}$ under tight assumptions. Their theorem translates directly into Theorem \ref{thm:fclt} for finite i.i.d s.b.p as discussed in Section \ref{subsec:fclt}.
Various versions of results in Section \ref{sec:asymp}, including Theorem \ref{thm:fclt} are also known in the successive sampling community \cite{holst, rosen, gordon83, bnw, sen}. For example, Bickel, Nair and Wang \cite{bnw} proved Theorem \ref{thm:gc} with convergence in probability when $E_n$ and $F$ have the same discrete support on finitely many values.

\section{Poisson coupling of size-biased permutation and order statistics}\label{sec:ppp}

Comparisons between the distribution of induced order statistics and order statistics of the same sequence have been studied in the literature \cite{dayang77, yang77, nagada94, hall94}. However, finite i.i.d s.b.p has the special feature that there exists an explicit coupling between these two sequences as described in Proposition ~\ref{prop:coupling}. Using this fact, we now derive Theorem \ref{thm}, which gives a Poisson coupling between the \emph{last} $k$ size-biased terms $X_n^{rev}[1], \ldots, X_n^{rev}[k]$ and the $k$ smallest order statistics $X_n^{\uparrow}(1), \ldots X^{\uparrow}_n(k)$ as $n \to \infty$. \footnote{The superscript `rev' indicates that the order statistics in consideration are arranged in increasing order, as opposed to decreasing, which has been the convention of this paper.}
The existence of a Poisson coupling should not be surprising, since it is well-known that the increasing sequence of order statistics $(X_n^{\uparrow}(1), X_n^{\uparrow}(2), \ldots)$ converges to points in a Poisson point process (p.p.p) whose intensity measure depends on the behavior of $F$ near the infimum of its support, which is $0$ in our case. This standard result in order statistics and extreme value theory dates back to R\'{e}nyi \cite{renyi}, and can be found in \cite{evtHaan}.

Our theorem is closely related to a result in PPY \cite[\S 4]{pitmanyor}. When $(X^\downarrow(1), X^\downarrow(2), \ldots)$ are ranked jumps of a subordinator, these authors noted that one can couple the size-biased permutation with the order statistics via the following p.p.p
\begin{equation}\label{eqn:ppois}
 N(\cdot) := \sum_{k\geq 1}\mathbf{1}[(X[k],Y^{\uparrow}(k)) \in \cdot] = \sum_{k\geq1}\mathbf{1}[(X^{\downarrow}(k),Y_k) \in \cdot]. \end{equation}
 $N(\cdot)$ has measure $ m(dxdy) = xe^{-xy}\, \Lambda(dx) dy. $ The first expression in (\ref{eqn:ppois}) defines a scatter of $(x,y)$ values in the plane listed in increasing $y$ values, and the second represents the same scatter listed in decreasing $x$ values. Since $\sum_{i\geq 1}X^\downarrow(i)$ is a.s. finite, the $x$-marginal of the points in  (\ref{eqn:ppois}) has the distribution of the size-biased permutation of the infinite sequence $(X^\downarrow(1), X^\downarrow(2), \ldots)$ since it prescribes the joint distribution of the first $k$ terms $X[1], \ldots, X[k]$ for any finite $k$.
PPY use this p.p.p representation to generalize size-biased permutation to $h$-biased permutation, where the `size' of a point $x$ is replaced by an arbitrary strictly positive function $h(x)$; see \cite[\S 4]{pitmanyor}.

More generally, one obtains a random permutation of $\mathbb{N}$ from ranking points $(x_i, y_i)$ in a Poisson scatter $N(\cdot)$ on $(0, \infty)^2$ according to either their $x$ or $y$ values. For $i = 1, 2, \ldots$, let $(x^\ast_i)$ and $(y^\ast_i)$ denote the induced order statistics of the sequence $(x_i)$ and $(y_i)$ obtained by ranking points by their $y$ and $x$ values in increasing order, respectively. For $j, k = 1, 2, \ldots$, define sequences of integers $(K_j), (J_k)$ such that $x^{\uparrow}(J_k) = x^{\ast}_k$, $y^{\uparrow}(K_j) = y^{\ast}_j$; see Figure \ref{fig:scatter}.

\begin{figure}[hbtp]
\vskip-1.5em
\includegraphics[width = 0.7\textwidth]{./pic1.pdf}
\vskip-4.5em
\caption{A point scatter on the plane. Here $J_1 = 5, J_2 = 3, J_3 = 1, J_4 = 4, J_5 = 2$, and $K_1 = 3, K_2 = 5, K_3 = 2, K_4 = 4, K_5 = 1$. The permutations $J$ and $K$ are inverses. Conditioned on $x^\uparrow(2) = a$ and $y_2^\ast = b$, the number of points lying in the shaded region determines $J_2 - 1$.}
\label{fig:scatter}
\end{figure}
Suppose $N(\cdot)$ has intensity measure $m$ such that for all $x, y \in (0, \infty)$, $m((0, x) \times (0, \infty)) < \infty, m((0, \infty) \times (0, y)) < \infty$. For $j \geq 1$, conditioned on $x(j) = x, y^{\ast}_j = y$,
\begin{equation}\label{eqn:kj}
K_j - 1 \stackrel{d}{=} Poisson\left(m((x, \infty) \times (0, y))\right) + Binomial\left(j-1, \frac{m((0, x) \times (0, y))}{m((0, x) \times (0, \infty))}\right),
\end{equation}
where the two random variables involved are independent. Similarly, for $k \geq 1$, conditioned on $x^\ast_k = x$, $y(k ) = y$,
\begin{equation}\label{eqn:jk}
J_k - 1 \stackrel{d}{=} Poisson\left(m((0,x) \times (y,\infty))\right) + Binomial\left(k-1, \frac{m((0, x) \times (0, y))}{m((0, \infty) \times (y, \infty))}\right),
\end{equation}
where the two random variables involved are independent. When $m$ is a product measure, as is the case of i.i.d s.b.p, it is possible to compute the marginal distribution of $K_j$ and $J_k$ explicitly for given $j, k \geq 1$. We demonstrate such computations in Proposition ~\ref{prop:compute.general}.

Before stating the theorem we need some technical results. The distribution of the last few size-biased picks depends on the behavior of $F$ near $0$, the infimum of its support. We shall consider the case where $F$ has `power law' near $0$, like that of a Gamma distribution.
\begin{lem} \label{lem:approx.G.by.Gamma}
Suppose $F$ is supported on $(0, \infty)$ with Laplace transform $\phi$. Let $u = \phi(y)$, $X_u$ a random variable distributed as $G_u(dx)$ defined in (\ref{eqn:Gu}). For $\lambda, a > 0$,
\begin{equation}
F(x) \sim \frac{\lambda^a x^a}{\Gamma(a + 1)} \mbox{ as } x \to 0, \label{eqn:assump.F}
\end{equation}
if and only if,
\begin{align}
\phi(y) & \sim \lambda^a/y^{a} \hspace{1em} \mbox{ as } y \to \infty \label{eqn:phi.asymp}.
\end{align}
Furthermore, (\ref{eqn:assump.F}) implies
\begin{align}
u^{-1/a}X_u & \stackrel{d}{\longrightarrow}  gamma(a+1,\lambda) \hspace{1em} \mbox{ as } u \to 0.
\label{eqn:xu.gamma}
\end{align}

\end{lem}
\begin{proof} The equivalence of (\ref{eqn:assump.F}) and (\ref{eqn:phi.asymp}) follows from a version of Karamata Tauberian Theorem \cite[\S 1.7]{regvar}. Assume (\ref{eqn:assump.F}) and (\ref{eqn:phi.asymp}). We shall prove (\ref{eqn:xu.gamma}) by looking at the Laplace transform of the non-size-biased version $X_u'$, which has distribution $F_u$. For $\theta \geq 0$,

\begin{equation}\label{eqn:laplace.xpu}
 E(\exp(-\theta X_u')) = \int_0^\infty u^{-1}\exp(-yx - \theta x)\, F(dx) = u^{-1}\phi(y + \theta) = \frac{\phi(y+\theta)}{\phi(y)}.
 \end{equation}
Now as $y \to \infty$ and $u = \phi(y) \to 0$, for each fixed $\eta > 0$, (\ref{eqn:phi.asymp}) implies
$$
E(\exp(-\eta u^{-1/a}X_u')) = \frac{\phi(y + \eta\phi(y)^{-1/a})}{\phi(y)}
\sim \frac{\lambda^a(y+\eta \lambda^{-1} y)^{-a} }{\lambda^ay^{-a}} = \left(\frac{\lambda}{\lambda+\eta}\right)^a.
$$
That is to say
\begin{equation}\label{eqn:xuprime.gamma}
u^{-1/a} X_u' \stackrel{d}{\longrightarrow} gamma(a, \lambda).
\end{equation}
Since $\phi$ is differentiable, (\ref{eqn:laplace.xpu}) implies $E(X'_u) = \phi'(y) / \phi(y)$. Now $\phi$ has an increasing derivative $\phi'$, thus (\ref{eqn:phi.asymp}) implies $\phi'(y) \sim a\lambda^a/y^{a+1}$ as $y \to \infty.$ Therefore,
$$u^{-1/a}E(X_u') = \frac{\phi'(y)}{\phi(y)^{1 + 1/a}} \to \frac{a}{\lambda},$$
which is the mean of a $gamma(a,\lambda)$ random variable. Thus the random variables $u^{-1/a}X'_u$ are uniformly integrable, so for any bounded continuous function $h$, we can compute
$$E(h(u^{-1/a}X_u)) = \frac{E[(u^{-1/a}X_u')h(u^{-1/a}X_u')]}{u^{-1/a}E(X_u')} \to \frac{E[\gamma_{a,\lambda}h(\gamma_{a,\lambda})]}{E(\gamma_{a,\lambda})} = E(h(\gamma_{a+1,\lambda})) $$
where $\gamma_{b,\lambda}$ is a $gamma(b,\lambda)$ random variable. This proves (\ref{eqn:xu.gamma}).
\end{proof}

Now suppose $F$ satisfies (\ref{eqn:assump.F}) for some $\lambda, a > 0$. By standard results in order statistics \cite[Theorem 2.1.1]{evtHaan}, as $n \to \infty$,
\begin{equation}\label{eqn:fidi.order}
\left(n^{1/a} X^{rev}_n(k), 1 \leq k \leq n\right) \stackrel{f.d.d.}{\longrightarrow} \left(\xi^{rev}(k): 1 \leq k < \infty\right),
\end{equation}
where f.d.d. stands for finite dimensional distribution. Here
\begin{equation}\label{eqn:fidi.sbs}
\xi^{rev}(k) = S_{(k)}^{1/a}\Gamma(a+1)^{1/a}/\lambda
\end{equation}
for $S_{(k)} = \epsilon_1 + \ldots + \epsilon_k$, where $\epsilon_i$ are i.i.d standard exponentials.
We now present the analogue of (\ref{eqn:fidi.sbs}) for the last few size-biased picks $X_n^{rev}[1], \ldots, X_n^{rev}[k]$ and the promised Poisson coupling.
\begin{thm} \label{thm}
Suppose that (\ref{eqn:assump.F}) holds for some $\lambda, a > 0$. \\
\textbf{a)} As $n \to \infty$,
\begin{equation}\label{15.a}
\left(n^{1/a} X_n^{rev}[k], 1 \leq k \leq n\right) \stackrel{f.d.d}{\longrightarrow} \left(\xi^{rev}[k]: 1 \leq k < \infty\right),
\end{equation}
with
\begin{equation} \label{15.b}
\xi^{rev}[k] = T_{(k)}^{1/a}\gamma_k / \lambda,
\end{equation}
 where $T_{(k)} = \epsilon'_1 + \ldots + \epsilon'_k$ for i.i.d standard exponential r.v. $\epsilon'_i$, $\gamma_k, k = 1, \ldots, n$ are i.i.d $gamma(a+1,1)$, and the $(\epsilon_k')$ and $(\gamma_k)$ are independent.
\vskip12pt \noindent
\textbf{b)} Let $\xi^{rev}(1) < \xi^{rev}(2) < \ldots $ be the increasing order statistics of the $\xi^{rev}[k]$'s. Then the $\xi^{rev}(k)$'s are jointly distributed as (\ref{eqn:fidi.sbs}), for a sequence $(\epsilon_k)$ which is implicitly defined in terms of the $(\epsilon'_k)$ and $(\gamma_k)$ above, and the f.d.d. convergence in (\ref{eqn:fidi.order}) and (\ref{15.a}) hold jointly.
\vskip12pt \noindent
\textbf{c)} For each $n$, let $J_n = (J_{nk}, 1 \leq k \leq n)$ be the permutation of $\{1, \ldots, n\}$ defined by
$$ X^{rev}_n[k] = X_n(J_{nk}). $$
As $n \to \infty$,
\begin{equation}\label{15.g}
\left(J_{nk}, 1 \leq k \leq n\right) \stackrel{f.d.d.}{\longrightarrow} \left(J_k: 1 \leq k < \infty\right),
\end{equation}
where $J_k$ is the random permutation of $\{1, 2, \ldots, \}$, defined by
\begin{equation}\label{15.h}
 \xi^{rev}[k] = \xi^{rev}(J_k)
 \end{equation}
 for $k = 1, 2, \ldots,$ and the f.d.d. convergence in (\ref{eqn:fidi.order}), (\ref{15.a}), (\ref{15.g}) all hold jointly.
\end{thm}
\begin{proof} By Lemma \ref{lem:approx.G.by.Gamma}, it is sufficient to prove the theorem for the case $F$ is $gamma(a, \lambda)$. Part a follows immediately from Proposition \ref{prop:newgb} and law of large numbers, in the same way that part b follows from the proof of (\ref{eqn:fidi.sbs}). For the last part, note that in principle everything has been presented as a function of the variables $T_{(k)}$ and $\gamma_k$. For $a, \lambda > 0$, define
$$\Psi_{a,\lambda}(s) = s^{1/a}\Gamma(a+1)^{1/a}/\lambda.$$
Observe that $\Psi_{a,\lambda}$ has inverse function $\Psi_{a,\lambda}^{-1}(x) = \lambda^a x^a / \Gamma(a+1)$. Rewrite (\ref{eqn:fidi.sbs}) as
\begin{equation}\label{eqn:xik}
\xi^{rev}(k) = \Psi_{a,\lambda}(S_{(k)}),
\end{equation}
then (\ref{15.b}) implies
\begin{equation}\label{16.a}
 S_k = \Psi_{a,\lambda}^{-1}(\xi^{rev}[k]) = \lambda^a \xi^{rev}[k]^a/\Gamma(a+1) = T_{(k)}\gamma_k^a / \Gamma(a+1).
 \end{equation}
Comparing (\ref{eqn:xik}) and (\ref{16.a}) gives a pairing between $S_k$ and $S_{(k)}$ via $\xi(k)$ and $\xi^{rev}[k]$. Hence we obtain another definition of $J_k$ equivalent to (\ref{15.h}):
$$ S_k = S_{(J_k)}. $$
Let $T_j$ be the $T$ value corresponding to the order statistic $S_{(j)}$ of the sequence $S_k$. That is,
$$T_j = T_{(K_j)}$$
where $(K_j)$ is a random permutation of the positive integers. By (\ref{16.a}), $(J_k)$ is the inverse of $(K_j)$. Together with (\ref{eqn:fidi.order}) and (\ref{15.a}), this implies (\ref{15.g}), completing the proof of part c.
\end{proof}
 Since the $T_{(k)}$'s are increasing ordered points of a p.p.p on $(0, \infty)$ with rate 1 and independent of the $\gamma_k$'s, the point process on the positive quadrant $(0, \infty)^2$
\begin{equation}\label{16.b}
 N(\cdot) = \sum_k\mathbf{1}[(S_k, T_{(k)}) \in \cdot] = \sum_j\mathbf{1}[(S_{(j)}, T_j) \in \cdot]
 \end{equation}
is a p.p.p for the measure $\mu$ with density
\begin{equation}\label{eqn:measure.mu}
 \frac{\mu(ds\, dt)}{ds\, dt} = \frac{1}{a} \Gamma(a+1)^{1/a}(s/t)^{1/a}\exp\{-(\Gamma(a+1)s/t)^{1/a}\}
 \end{equation}
The random permutation $J_k$ and its inverse $K_j$ are obtained from the p.p.p in (\ref{16.b}) in the $(s,t)$ plane by ranking the points in two different orders according to $s$-value or $t$-value. Since the projection of a Poisson process is Poisson, the $s$ and $t$-marginal of $\mu$ is just Lebesgue measure. Thus, the ranked $s$-values form a p.p.p with rate 1 which determine the ranked $\xi(j)$ by the \emph{deterministic} increasing transformation
$$ \xi^{rev}(k) = \Psi_{a,\lambda}(S_{(k)}).$$
Furthermore,
$$\xi^{rev}[k] = \xi^{rev}(J_k) = \Psi_{a,\lambda}(S_k)$$
comes from the same set of $s$-values listed in order of increasing $t$-values.

Marginal distributions of the random permutation $(J_k)$ and its inverse $(K_j)$ are given in (\ref{eqn:kj}) and (\ref{eqn:jk}). Note that for $k = 1,2, \ldots$, $S_k = T_{(k)}\gamma_k^a/\Gamma(a+1)$ for i.i.d $\gamma_k$ distributed as $gamma(a+1,1)$, independent of the sequence $(T_{(k)})$, and $T_k = \Gamma(a+1)S_{(k)}\tilde{\epsilon}_k^{-a}$ for i.i.d standard exponentials $\epsilon_k$, independent of the sequence $(S_{(k)})$ but not of the $\gamma_k$'s. Thus by conditioning on either $S_{(k)}$ or $T_{(k)}$, one can evaluate (\ref{eqn:kj}) and (\ref{eqn:jk}) explicitly. In particular, by a change of variable $r = \Gamma(a+1)^{1/a}(s/t)^{1/a}$, one can write $\mu$ in product form. This leads to the following.
\begin{prop}\label{prop:compute.general}
For $j \geq 1$, conditioned on $S_{(j)} = s, T_j = \Gamma(a+1)\, s \, r^{-a}$ for some $r > 0$, $K_j-1$ is distributed as
\begin{equation}\label{eqn:marginal.kj}
Poisson(m(s,r)) + Binomial(j-1, p(s,r))
\end{equation}
with
\begin{equation}\label{eqn:ppp.density.kj}
m(s,r) = asr^{-a}\, \int_r^\infty x^{a-1}e^{-x}\, dx,
\end{equation}
and
\begin{equation}\label{eqn:psr}
p(s,r) = as^{2/a-2}r^{a-2}\int_0^r x^{a-1}e^{-x}\, dx,
\end{equation}
where the Poisson and Binomial random variables are independent.
Similarly, for $k \geq 1$, conditioned on $T_{(k)} = t, S_k = t \, r^a / \Gamma(a+1)$ for some $r > 0$, $J_k - 1$ is distributed as
\begin{equation}\label{eqn:marginal.jk}
 Poisson(m'(t, r)) + Binomial(k-1, p'(t,r))
 \end{equation}
with
\begin{equation}\label{eqn:ppp.density.jk}
m'(t,r) = t \, \left(\frac{r^a + a \, \int_r^\infty x^{a-1}e^{-x}\, dx}{\Gamma(a+1)} - 1\right),
\end{equation}
and
\begin{equation}\label{eqn:ptr}
p'(t,r) = \Gamma(a+1)^{1-2/a}\, at^{2/a-2}r^{a-1/a}\int_0^r x^{a-1}e^{-x}\, dx,
\end{equation}
where the Poisson and Binomial random variables are independent.
\end{prop}
\begin{prop}[Marginal distributions of $K_1$ and $J_1$]\label{prop:compute}
Suppose that (\ref{eqn:assump.F}) holds for some $\lambda > 0$ and $a = 1$. Then the distribution of $K_1$, the $k$ such that $\xi(1) = \xi^{rev}[k]$, is a mixture of geometric distributions, and so is that for $J_1$, the $j$ such that $\xi^{rev}[1] = \xi(j)$. In particular,
\begin{equation}\label{18.k1}
P(K_1 = k) = \int_0^\infty p_rq_r^{k-1}e^{-r}\, dr
\end{equation}
where $p_r = r/(r + e^{-r})$, $q_r = 1 - p_r$, and
\begin{equation}\label{18.j1}
P(J_1 = j) = \int_0^\infty \tilde{p}_r\tilde{q}_r^{j-1}re^{-r}\, dr
\end{equation}
where $\tilde{p}_r = 1/(r+e^{-r})$, $\tilde{q}_r = 1 - \tilde{p}_r$.
\end{prop}
\begin{proof} When $a = 1$, $\int_r^\infty t^{a-1}e^{-t}\, dt = e^{-r}$. Substitute to (\ref{eqn:ppp.density.kj}) and (\ref{eqn:ppp.density.jk}) give
$$m(s,r) = sr^{-1}e^{-r}, \hspace{1em} m'(t,r) = t(r-1+e^{-r}).$$
By a change of variable, (\ref{eqn:measure.mu}) becomes
$$\frac{\mu(ds\, dr)}{ds\, dr} = se^{-r}, \hspace{1em} \frac{\mu(dt\, dr)}{dt\, dr} = tre^{-r}.$$
Thus, \emph{conditioned on $s$ and $r$}, $K_1 - 1$ is distributed as the number of points in a p.p.p  with rate $r^{-1}e^{-r}$ before the first point in a p.p.p with rate 1. This is the geometric distributions on $(0, 1, \ldots)$ with parameter $p_r = 1/(1 + r^{-1}e^{-r})$. Since the marginal density of $r$ is $e^{-r}$, integrating out $r$ gives (\ref{18.k1}). The computation for the distribution of $J_1$ is similar.
\end{proof}
One can check that (\ref{18.k1}) and (\ref{18.j1}) sum to 1. We conclude with a `fun' computation. Suppose that (\ref{eqn:assump.F}) holds for some $\lambda > 0$ and $a = 1$. That is, $F$ behaves like an exponential c.d.f near $0$. By Proposition \ref{prop:compute}, $E(J_1) = 9/4$ and $E(K_1) = \infty$. That is, the last size-biased pick is expected to be almost the second smallest order statistic, while the smallest order statistic is expected to be picked infinitely earlier on in a successive sampling scheme(!). The probability that the last species to be picked in a successive sampling scheme is also the one of smallest species size is
\begin{align*}
\lim_{n \to \infty} P(X_n^{rev}[1] = X_n(1)) &= P(\xi^{rev}[1] = \xi(1)) = P(J_1 = 1) = P(K_1 = 1) \\
&= \int_0^\infty \frac{re^{-r}}{r + e^{-r}}\, dr = 1 - \int_0^1 \frac{u}{u - \log u}\, du \approx 0.555229
\end{align*}

\bibliographystyle{plain}