\section{Introduction}

Complex optical systems with several optical elements are challenging to cool down while keeping aligned due to thermal contraction and changes in refractive index \cite{Purdy2012}. Therefore such systems are often based on optical free space access cryostats \cite{Purdy2012, Kuhn2014, Peterson2016, Rossi2018} to allow for external compensation of the optical misalignment. Alternatively, fiber coupled optical systems with cryogenic compatible actuators can be used leading to an increase in design complexity and requiring active control protocols while cooling down \cite{Verhagen2012, Fogliano2021}. Such systems are of particular interest for quantum optomechanical experiments that require low temperature environment and high collection efficiencies \cite{Aspelmeyer2014, Cohen2015, Riedinger2016, Galinskiy2020, Fedoseev2021} because the rate of thermal decoherence is proportional to the temperature of the environment \cite{Aspelmeyer2014}. In quantum optomechanical experiments detection of a scattered photon projects the mechanical system into a specific quantum state. Therefore the total collection efficiency $\eta$ of the heralding photons is essential for high fidelity operations and affects the duration of data collection $\tau$ in such experiments:  $\tau\propto 1/\eta^n$, $n\geq2$ \cite{Riedinger2018, Weaver2018a}. A constituent of the total collection efficiency is the coupling efficiency of the cryogenic optical cavity to the outside optical path, often an optical fiber. In a Helium-4 flow cryostat with free space optical access a coupling efficiency of 60\% is reported \cite{Galinskiy2020}, while for a free space optical access dilution cryostat a coupling efficiency of 85\% was achieved \cite{Peterson2017}. Keeping mode matching between a single-mode fiber and a free space cavity with cryogenic actuators is even more challenging, requiring many degrees of freedom to be controlled during cooling down. The resulting mode matching in such systems is quite limited, with reported values in the range of 30-40\% \cite{Buters2017, Xu2019}.

Here we present a fiber coupled free space 10 cm long Fabry-Perot resonator which is designed to be fully axially symmetric. This allows for room temperature prealignment of the system such that the optimal alignment is obtained at cryogenic temperatures. Without any position adjustments during a cool down 96$\pm$2\% coupling efficiency between a single-mode fiber and the optical cavity with a finesse of 12000 at the optical wavelength of 1064 nm is achieved. The elimination of actuators  does not only simplify the system but also significantly improves the mechanical stability which is of vital importance for optomechanical experiments \cite{Purdy2012}.

\section{Setup}

Figure 1 shows the design concept of our custom-built optical setup.  The main body and all the elements in contact with the optical components are made of Invar

\begin{figure}[ht]
\centering
\includegraphics[scale=0.37]{fig1.eps}
\caption{Cryogenic setup. (a) The setup is axially symmetric and is made mostly of Invar 36. The fiber ferrule is held in a differential screw, the part shown in red can be rotated, the inner part of the screw is pushed with a leaf spring and allows for purely translational shifts of the ferrule along the optical axis. The lenses and mirrors are held in compartments (light grey) which are pushed against the setup modules (dark grey) with springs. These compartments can be shifted perpendicularly to the optical axis by 4 screws (right sketch). The screws are retracted after the alignment and the compartments are held by friction during cooling down. (b) The setup is placed in a dilution cryostat. The fiber tip of the cryogenic setup is not angled and is not antireflection-coated, the light reflected from the tip interferes with the light both reflected and leaking from the cavity.}
\end{figure}

\begin{table*}[ht]
\caption{\label{tab:table3} Positions and properties of the optical elements calculated using Gaussian beam propagation.}
\begin{ruledtabular}
\begin{tabular}{m{5cm}cccc}
    & lens1 & lens2 & mirror1 & mirror2 \\
 \hline
 Focal length, mm & 4.5 & 48 & 25 & 25 \\
 Distance from the fiber tip, mm & 4.2 & 40 & 94 & 190 \\
 Material & D-ZK3 & N-BK7 & fused silica & fused silica\\
\end{tabular}
\end{ruledtabular}
\end{table*}

\begin{figure*}[!ht]\centering
\includegraphics[scale=0.8]{fig2.eps}
\caption{Total reflectivity $|r_\mathrm{total}|^2$ far away from cavity resonances. (a) Theoretical interference of the light back-reflected from the fiber tip and mirror1 coupled into the fiber vs fiber coupling efficiency of the light back-reflected from mirror1 $\alpha_m$. The red (blue) line is the maximum (minimum) total reflectivity $|r_\mathrm{total}|^2$ for a given coupling efficiency. (b) Red (blue) lines/circles show the measured maximum (minimum) reflectivity due to the interference described in (a) during cooling down of the setup. The lines show the values extracted from measured visibility of the interference, while the circles show the reflectivity calculated as the ratio of light intensities measured by the photodetector and the powermeter with the ratio scaling corresponding to the room temperature measurements. (c) Coupling efficiency of back-reflected light from mirror1 $\alpha_m$ vs temperature extracted from the measured interference visibility. In the inset we demonstrate the reversibility of the coupling efficiency changes by showing the coupling efficiency $\alpha_m$ vs temperature for cooling down and warming up. The data is from another experimental run with a different optics positioning, magenta circles - cooling down, orange circles - warming up.
}
\end{figure*}

\noindent 36 to minimize thermal contraction. The light reaches the setup via an optical fiber (1060XP, numerical aperture 0.14) which is terminated with a ceramic ferrule, not angled. The ferrule is held by a differential screw (Invar 36) with a pitch difference of 0.1 mm, where the rotatable part is depicted in red. The inner part of the differential screw is connected to a leaf spring (Invar 36), which allows translation along the optical axis without rotation, and eliminates movement of the inner part due to the tolerances in the threading. Tests of the differential screw showed that the light originating from the fiber ferrule and passing through the first lens only deviates by no more than 1 mm at a distance of 1 m from the ferrule when the screw is rotated by one turn.

Another important design feature is that the lenses and mirrors are held in compartments (light grey) which are pressed towards the body of the setup (dark grey) along the optical axis with springs. The compartments can be shifted perpendicularly to the optical axis by four screws (80 PTI). When the alignment is finished, the screws are retracted and removed which does not change the position of the compartments, and the compartments are held in place by friction. A small shift of the lenses perpendicularly to the optical axis leads to small changes of the position and direction of the light between lens2 and mirror1, effectively acting as a periscope.

The system is designed to be fully axially symmetric: the thermal contraction and change of optical properties of the elements lead to the shift of the waist position of the cavity mode relative to the fiber tip only along the mechanical axis, which can be precompensated before cooling down. The system can be aligned in such a way that the light in the setup travels along the mechanical axis of symmetry. In a conventional design with fixed mirrors and a periscope the cavity mirrors are, to some degree, randomly positioned not guaranteeing the coincidence of the optical and mechanical axes. This issue is especially important for a cavity close to concentric (as here) where the angle between the two axes in the cavity $\alpha$ is given by $\tan{\alpha}=e/(2R-L)$ in the symmetric cavity case, $e$ is the error of positioning of the mirrors, $R$ is the mirrors radius of curvature, $L$ is the distance between the mirrors. For our case of $L=94$ mm, $R=50$ mm and a typical error of $e=0.1$ mm would lead to $\alpha=0.017$ rad or 0.8 mm shift of the beam center on the mirror.

\section{Alignment procedure}
Positions of the optical elements leading to the theoretically perfect mode matching are shown in Table I. The aligning procedure starts with positioning mirror1 on the mechanical axis of the setup (z-axis). Light going from the fiber through lens1 and a pinhole at the position of mirror2 on the mechanical axis is imaged with an intensity calibrated camera. This allows to position lens1 such that the optical and mechanical axes coincide. Next mirror1 is added and its position is adjusted to maximize the light intensity coupled back to the fiber. This defines the final position of mirror1.

Next, lens2 is placed at the calculated optimal position (see Table I). The fiber ferrule z-position, (x,y) positions of lens1 and lens2 are adjusted to maximize the light intensity coupled back to the fiber. The fiber tip is not anti-reflection coated and is not angled, and 3.4\% of light is reflected back into the fiber. This reflection interferes with the light reflected from mirror1. The visibility of this interference allows to calculate the coupling efficiency of the light back-reflected from mirror1 [see Fig. 2(a)]. The fine adjustment of the fiber ferrule z-position and the lens positions (ferrule z-position, (x,y) positions of lens1 and lens2) are made to minimize the interference visibility.

Now, mirror2 is added. First, its position is changed to maximize the transmission of the 0,0 transverse mode, fine adjustment is achieved by minimizing the 0,1 and 1,0 transverse Hermite-Gaussian modes transmission. The difference in the beam diameters of the light coupled to the cavity and the cavity mode is best visible in the non-zero transmission through 0,2 and 2,0 transverse Hermite-Gaussian modes. This gives an indication whether lens2 is positioned properly.

Lens2 can be shifted along the z-axis and is fixed by two threaded retaining rings. Final position of lens2 is found through repeating the two previous steps: for a set of different z positions of lens2 z-position of the ferrule and (x,y) positions of lens1 and lens2 are adjusted to maximize the light reflected from mirror1 and coupled to the fiber, and mirror2 position is adjusted to minimize the 0,1 and 1,0 couplings. The lens2 z-position with minimum 0,2 and 2,0 transmission is the optimal one. In the end of this procedure we achieved a configuration where 98.5\% of light is transmitted through the 0,0 mode and the rest is almost equally distributed between modes 0,1, 1,0, 0,2 and 2,0 modes. Transmission through other higher-order modes is at least one order of magnitude lower.

Finally, the setup is precompensated for the changes due to cooling down: in the ideal axially symmetric case the change of the optical properties of the components and thermal contraction effectively result in a shift of the waist z-position of the cavity mode at the fiber and its spot size. The spot size change has negligible effect on the fiber-cavity mode overlap. This makes it possible to precompensate for the cooling down effects only by shifting the fiber ferrule along the optical axis outwards from the cavity relative to the optimal room temperature coupling configuration. (x,y) positions of the lenses are readjusted to maximize the coupling of the light reflected from mirror1 to the fiber after the ferrule shift to compensate for small beam-pointing changes mentioned above. In the experiment reported here, the ferrule was retracted by 32$\pm$3 $\mu$m before the cool down. The retraction for this cool down is an estimate of the required precompensation based on the two previous cool downs. Each cool down provides a better estimate of the required precompensation for the next cool down.

\section{Measurements}
The fiber tip is not anti-reflection coated and is not angled, therefore the reflection signal measured by the photodetector is the result of interference of the reflections from the fiber tip, mirror1 and mirror2. This allows us to accurately measure the change of the coupling efficiency of the light back-reflected from mirror1 to fiber mode $\alpha_m$, and coupling efficiency of the cavity mode to fiber mode $\alpha_c$. In this article, we define "coupling efficiency" of a field $E$ as the intensity mode overlap of the fiber field $E_f$ and $E$ diminished by the losses $\epsilon$:

\begin{equation}
\alpha=\frac{|\int E^*E_fdA|^2}{\int E^*_fE_fdA\int E^*EdA}(1-\epsilon)
\end{equation}

In particular, if the total intensity losses on the way from the fiber tip to mirror1 are $\epsilon_\mathrm{single}$ (say, a dust particle on the fiber tip and not mode-matched reflections from lens1, lens2 and back mirror1 surfaces) then

\begin{equation}
\begin{aligned}
\alpha_m&=\frac{|\int E_m^*E_fdA|^2}{\int E^*_fE_fdA\int E^*_mE_mdA}(1-\epsilon_\mathrm{single})^2\\
\alpha_c&=\frac{|\int E_c^*E_fdA|^2}{\int E^*_fE_fdA\int E_c^*E_cdA}(1-\epsilon_\mathrm{single})
\end{aligned}
\end{equation}

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.8]{fig3.eps}
\caption{Total reflectivity $|r_\mathrm{total}|^2$. (a) Room temperature light frequency scan in the case of optimal alignment before the precompensational shift of the ferrule. The cavity resonance dips were measured separately. (b) A particular cavity resonance is traced during cooling down (after the precompensation). The black line shows the off-resonant reflectivity in the vicinity of the traced resonance, red (blue) - maximum (minimum, excluding cavity resonances) reflectivity when the laser frequency is swept over 1 FSR of the fiber tip - mirror1 system. The insets show the relative phases between the fiber tip - mirror1 system and the cavity. Green - resonance reflectivity of the traced cavity resonance. Black circles show the resonance reflectivity values mapped from the room temperature alignment shown in a). The very similar reflection properties of the optimal alignment at room temperature and at 130 K shows that the setup alignment at 130 K is very close to the best achievable alignment at room temperature. (c) Black, red, blue and green - same as in b). Black squares represent theoretical reflectivity dip values calculated with $\alpha_c=94$\%.
}
\end{figure*}
The expression for $\alpha_m$ includes the term $(1-\epsilon_\mathrm{single})^2$ as the light traveling from the fiber and back is experiencing the losses $\epsilon_\mathrm{single}$ twice. The total reflectivity amplitude is given by

\begin{equation}
\begin{aligned}
&r_\mathrm{total}=-r_1+\frac{(1-r_1^2)Me^{i\phi_1}}{1-r_1Me^{i\phi_1}} \mathrm{,  where}\\
&M=\sqrt{\alpha_m}r_2-\frac{\alpha_c t_2^2r_3e^{i\phi_2}}{1-r_2r_3e^{i\phi_2}}
\end{aligned}
\label{eq_main}
\end{equation}

and $r_1, r_2, r_3$ are reflectivity amplitudes of the fiber tip, mirror1, mirror2, respectively. $t_2$ is transmissivity amplitude of mirror1, $\phi_1=2\pi/\lambda\times2L_\mathrm{t1}$ with $\lambda=1064$ nm and $L_\mathrm{t1}$ being fiber tip - mirror1 optical distance, $\phi_2=2\pi/\lambda\times2L_{c}$ with  $L_{c}$ being  cavity length. According to the specifications, mirror1 and mirror2 have an intensity transmissivity of 250 ppm; the scattering and absorption loss is not specified, but our measurements are consistent with 30 ppm. $r_1=\sqrt{0.034}$ based on the Fresnel coefficients for silica.

First, let's consider optical frequencies far away from the cavity resonances (detuning from the nearest cavity resonance $>$10 cavity linewidths, cavity finesse is 12000). For these frequencies the measured reflection signal is the interference of the system fiber tip - mirror1. In this case Eq. \ref{eq_main} can be simplified:

\begin{equation}
r_\mathrm{total}=\frac{-r_1+\sqrt{\alpha_m}r_2 e^{i\phi_1}}{1-r_1r_2 \sqrt{\alpha_m} e^{i\phi_1}}
\label{eq_simple}
\end{equation}

The theoretical total reflectivity $|r_\mathrm{total}|^2$ minimum and maximum values ($\phi_1=0$ and $\phi_1=\pi$ respectively) are shown in Fig. 2(a) as a function of coupling  efficiency of back-reflected light from mirror1 $\alpha_m$. For known $r_1$ and $r_2$ the visibility of this interference is a function of $\alpha_m$ only. This allows us to accurately estimate the coupling efficiency $\alpha_m$ by measuring the visibility of the interference between the fiber tip and mirror1 reflections during cooling down of the setup, no calibration is needed. The minimum and maximum of the total reflectivity can be calculated from the measured visibility data, these are plotted by the solid lines in Fig. 2(b).

The total reflectivity can also be calculated directly as a calibrated ratio of the reflection signal and the powermeter signal, see Fig. 1(b). Resulting minimum and maximum reflectivity values are shown by circles in Fig. 2(b). This estimation assumes that the transmission of the optical fiber going to the cryostat does not change with time/temperature. There is a fiber connector in the cryostat at the 4 K plate, which might explain the difference between the two estimates of the total reflectivity. The coupling efficiency of the light back-reflected from mirror1 $\alpha_m$ calculated using the visibility data is shown in Fig. 2(c). There is an initial discontinuity in $\alpha_m$ at room temperature. This change occurred during evacuation of the cryostat vacuum chamber causing a change in properties of the optical elements. $\alpha_m$ reaches 95\% at 140 K setting an upper bound on the total losses: $\epsilon_\mathrm{single}\leq0.025$. Lens1, lens2 and the back sides of the mirrors have anti-reflection coatings with a reflectivity of $\sim$0.2\% per surface (specs) setting the lower bound on the total losses: $\epsilon_\mathrm{single}\geq0.01$.

The change of the coupling efficiency with temperature is reversible: $\alpha_m$ is the same for a given temperature during cooling down and warming up as shown in the inset of Fig. 2(c) (data from previous cool down).

Next, we consider interference of the fiber tip, mirror1 and mirror2. The highest cavity mode overlap we could achieve at room temperature is 0.985 as stated in the section Alignment Procedure. A laser frequency scan covering 9 cavity resonances is shown in Fig. 3(a). A fit of Eq. \ref{eq_main} gives $\alpha_c=0.965$ for this data, providing an estimate of the losses $\epsilon_\mathrm{single}=0.02\pm0.005$ (0.985/0.965-1). The interference pattern of fiber tip - mirror1 results in the sinusoidal reflectivity variation. The optical path from the fiber tip to mirror1 is 6\% longer than the cavity length, which agrees well with the ratio of the cavity and fiber tip - mirror1 free spectral ranges of $8.5/8$ extracted from Fig. 3(a).

During the cooling down the laser frequency was actively controlled in a way to trace the same resonance of the cavity. To achieve this, the laser frequency was scanned over 100 cavity linewidths once a second and the laser frequency was corrected to keep the cavity resonance approximately in the middle of each frequency scan. The reflectivity on cavity resonance is shown by green lines in Fig. 3(b) and 3(c). Every 10 minutes the laser was scanned over 1 fiber tip - mirror1 FSR to obtain the minimum and maximum reflectivity (blue and red lines in Fig. 3(b) and (c)).

The data plotted by black, blue and red lines in Fig. 3(a) and (b) allows to estimate $\phi_1$ from Eq. \ref{eq_main}. This phase information allows to map the reflectivity dip values from the optimal room temperature alignment onto the $130\pm10$ K temperature interval as depicted by black circles in Fig. 3(b). To do this $\phi_1$ is estimated for each resonance from Fig. 3(a) and the resonance reflectivity values from Fig. 3(a) are plotted for temperatures corresponding to the same set of $\phi_1$ values estimated from Fig. 3(b). Minimum and maximum reflectivity values together with the reflectivity on cavity resonance are practically the same for the optimal alignment at room temperature (before precompensation) and at around 130 K which can be seen by direct comparison of the data from Fig. 3(a) and (b). Therefore, the cavity indeed can be precompensated to reach near optimal alignment at cryogenic temperatures.

During the two temperature intervals shown in Fig. 3(b) and (c) $\phi_1$ was changing by approximately $2\pi$. The data with $\phi_1$ spanning $2\pi$ is enough to extract a single value of $\alpha_c$ consistent with Eq. \ref{eq_main}: 0.96 for Fig. 3(b) and 0.94 for Fig. 3(c). This analysis was performed on 7 additional temperature intervals giving 7 additional single values $\alpha_c$, these estimates of the cavity coupling efficiency are shown by blue circles in Fig. 4.

\begin{figure
\centering
\includegraphics[scale=0.9]{fig4.eps}
\caption{Cavity mode coupling efficiency (left axis). Blue circles - fits to the measured reflection data. Red line - shift of the cavity mode waist at the fiber tip (see inset) calculated from the data of Fig. 2(c) under the assumption that the misalignment due to cooling down is axial symmetric. Such a shift would result in the cavity mode coupling efficiency depicted by the blue line. The arrow represents the ferrule shift made at room temperature to precompensate for the changes due to cooling down. Dashed lines are guides to the eye.
}
\end{figure}

Next, we check the hypothesis whether the cavity misalignment is axisymmetric during the cool down. Under this assumption the mode overlap of the back-reflected light from mirror1 with the fiber mode is a single variable function of the distance $\Delta z_m$ between the fiber tip and the waist of the mirror1 back-reflected light, depicted by the red arrows in the inset of Fig. 4. We estimated that the beam waist radius of the back-reflected light changes negligibly small during the cool down. The back-reflected light from mirror1 mode overlap with the fiber mode can be calculated as $\alpha_m/(1-\epsilon_\mathrm{single})^2$, see Eq. 2, and then used to calculate $\Delta z_m$ as all parameters of the Gaussian beams are known.

The distance between the fiber tip and the waist of the light leaked from the cavity at the fiber tip [see inset Fig. 4] $\Delta z_c$ is twice smaller than $\Delta z_m$ because the light leaked from the cavity goes through lens1 and lens2 once, while the light back-reflected from mirror1 goes through lens1 and lens2 twice: $\Delta z_c=\Delta z_m/2$. Calculated values of $\Delta z_c(T)$ are shown in Fig. 4 by the red line. This estimate $\Delta z_c(T=293\:\mathrm{ K})=30$ $\mu m$ agrees well with the precompensation shift of the fiber ferrule of 32$\pm$3 $\mu$m. Further, calculated values of $\Delta z_c(T)$ are used to calculate $\alpha_c(T)$ shown by the blue line in Fig. 4. These values are consistent with the previous estimates of $\alpha_c(T)$ based on the reflection dip measurements (blue circles) confirming the hypothesis of mainly axisymmetric misalignment of the setup during cooling down. The calculated fiber ferrule precompensation shift to reach the highest cavity mode coupling efficiency at mK temperatures is 45 $\mu$m.

Next, we discuss vibrational properties of the setup. The setup was cooled in a dilution cryostat with cooling provided by a pulse tube cryocooler known to introduce high levels of vibrations. We used a mass-spring system \cite{DeWit2019} acting as a fourth-order mechanical low pass filter with corner frequency of 50 Hz. The filter was placed between the 20 mK plate and the cryogenic cavity setup. The optical linewidth of the cavity was 130 kHz, the Pound-Drever-Hall lock was extremely stable with visible mechanical vibrations at 2-4 kHz just above the noise floor of our detection scheme while using 20 $\mu$W of probe laser intensity. Due to the monolithic design of the cryogenic cavity the most likely part susceptible to vibrations at 2-4 kHz are the mirror spring-loaded holders.

This setup is designed for quantum optomechanical experiments \cite{Fedoseev2021} where temperatures below 1 K and high collection efficiencies are required \cite{Riedinger2016, Galinskiy2020}. For this, we have made an opening in the center of the cavity where an Invar 36 based holder can be mounted with a Si chip with a transparent high-stress SiN membrane. The membrane can be positioned at the waist of the cavity mode (beam waist radius 60 $\mu$m) and vibrational modes of the membrane are coupled to the light inside the cavity. The cavity finesse with the membrane is $\sim$10000, requiring the membrane to maintain its orientation perpendicularly to the optical axis during cooling down within $\sim10^{-4}$ rad according to our experience. Cooling down the setup with the membrane shows the behavior of the reflectivity dip similar to Fig. 3(c), indicating that the loss due to membrane tilt is below $10^{-5}$ per membrane transmission. To maximize the collection efficiency of the photons scattered on the vibrations of the membrane, mirror2 will be exchanged for a 10 ppm transmission mirror to send the light leaking from the cavity predominantly towards the fiber. We would like to note that the absence of an anti-reflection coating on the fiber tip does not diminish the collection efficiency in the case of the ideal cavity - fiber mode matching and transmissivity of mirror1 being much larger than the transmissivity of mirror2: the scattered photon will leak from the cavity preferentially through mirror1 and will be either transmitted to the fiber or reflected from the fiber tip, and after entering the cavity will leak again towards the fiber.

\section{Conclusions}
We developed and characterized a single-mode fiber-coupled cryogenic optical cavity enabling cavity mode coupling efficiency above 90\% at mK temperatures. In particular, the setup does not require realignment during cooling down, due to the axisymmetric design and low expansion of the body made from Invar. The optics are aligned at room temperature to precompensate for the thermal contraction and change of refractive index during cooling down. The cavity mode coupling efficiency increases upon the start of cooling down from 81$\pm$2\%, reaches 96$\pm$2\% at 130 K and stabilizes at 92$\pm$2\% at dilution cryostat temperatures. As there are no actuators in the setup and due to its monolithic design, it has increased stiffness which reduces its susceptibility to vibrations in the cryostat.

The setup is designed to operate in the membrane-in-the-middle configuration \cite{Thompson2008} for optomechanical experiments in the quantum regime \cite{Fedoseev2021} where high collection efficiency is essential \cite{Galinskiy2020}.

\section{Acknowledgements}
We would like to thank K. Heeck for useful discussions of the practical design of the setup. This work was supported by the NSF Quantum Foundry through Q-AMASE-i program award number DMR-1906325 and NSF QuIC-TAQS OMA-2137740, and the NWO Gravitation-grant Quantum Software Consortium - 024.003.037.

\section{Author Declarations}
\subsection{Conflict of interest}
The authors have no conflicts to disclose.

\section{DATA AVAILABILITY}
The data that support the findings of this study are available from the corresponding author upon reasonable request.


\section{ Introduction }
For a fixed reflection group associated  with a root system $R$  and  for a multiplicity function $k\geq 0$, the   $(k,a)$-deformed  harmonic oscillator is given by
 \begin{equation*}
   \Delta_{k,a}= \|x\|^{2-a}\Delta_{k}-\|x\|^a
 \end{equation*}
  where  $a > 0$ is  a   parameter  and  $\Delta_k$ is   the Dunkl Laplacian operator on  $\mathbb{R}^d$.  This operator gives rise to  the  semigroup
 \begin{equation*}
   \mathscr{J}_{a}(z) = \exp\left(\frac{z}{a}\Delta_{k,a}\right)
 \end{equation*}
 for $z \in \mathbb{C}$ such that $ Re(z) \geq0$, first featured and studied in \cite{BKO}, where the authors defined in $L^2(\mathbb{R}^d,|x|^{a-2}\upsilon_{k,a}(x)dx )$ an unitary operator called  the $(k, a)$-generalized Fourier transform
  \begin{equation*}\label{}
     \mathscr{F}_{k,a}= e^{i\frac{\pi}{2}\frac{ d+a-2+\sum_{\alpha\in R} k(\alpha)}{a}}\mathscr{J}_{k,a}(\frac{i\pi}{2})
  \end{equation*}
which  can be    expressed as integral  transform:
  \begin{equation*}\label{ }
     \mathscr{F}_{k,a}(f)( \xi)= c_{k,a}\int_{\mathbb{R}^N} B_{k,a}(\xi,x) f(x)|x|^{a-2} \prod_{\alpha\in R}|\langle x,\alpha\rangle|^{k(\alpha)}dx.
  \end{equation*}
 with certain constant $c_{k,a}$. In particular, the case $a = 2$ corresponds to     Dunkl transform. Formal expressions   for $ B_{k,a}$    have been derived in \cite{BKO} as a  series representation,
but these expressions are not very useful from the analytic point of view.
 \par  in the one   dimensional case  the kernel $B_{k,a}$  is  given  by
 \begin{equation}\label{BB}
     B_{k,a}(\lambda,x)= \mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda x|^{a/2}\right)+
 m_{k,a}\lambda x \mathcal{J}_{\frac{2k+1}{a}}\left(\frac{2}{a}|\lambda x|^{a/2}\right), \qquad \lambda,x\in \mathbb{R}
 \end{equation}
 where
$$m_{k,a}=e^{\frac{-i\pi}{a}}\frac{\Gamma\left(\frac{2k+a-1}{a} \right)}{a^{\frac{ 2}{a}}\Gamma\left(\frac{2k+a+1}{a} \right)} $$
and $\mathcal{J}_{\nu}$  is the normalized  Bessel function.
\begin{eqnarray}\label{Jj}
\mathcal{J}_{\nu}(z)= \Gamma(\nu+1)\left(\frac{z}{2}\right)^{-\nu} J_\nu(z)=\Gamma(\nu+1)\sum_{n=0}^\infty \frac{(-1)^n\left(\frac{z}{2}\right)^{2n}}{ n!\Gamma(\nu+n+1)}.
\end{eqnarray}
Restricting then  to one dimensional case,  one of the classic problems that arises  is to describe  the product two $B_{k,a}'s$   in a most convenient way that is
$$B_{k,a}( \lambda,x)B_{k,a}(\lambda, y)=\int B_{k,a}(\lambda,z) d\gamma_{x,y}^{k,a}(z)$$
with $\gamma_{x,y}^{k,a}$ are  measures on $\mathbb{R}$
  which are uniformly
bounded with respect to total variation norm.
This  formula was established in  \cite{MSM}  for $a=\frac{2}{n}$,  $n\in \mathbb{N}$. The author's approach makes use of the   well-known  Gegenbauer's addition theorem for the Bessel functions. Our purpose
here is to  extend the formula of  \cite{MSM}  to the case $a>0$. To be more precise, $\gamma_{x,y}^{k,a}$ will
be derived  in terms of the associated Legendre functions which involved in the    infinite  integral of   product of three Bessel functions of the first kind, due to Macdonal \cite{Mac},(see also, \cite{Wat}). Through it,  and via  Hankel transform theory we present  some formulas for integrals involving  Bessel functions or their product.
 \section{Main Results }
 In this section, we establish two  integral formulas, which are expressed as Hankel transform of associate Legendre functions.
 \par    Recalling first
 the Macdonal  integral, that when $x$ and y are positive,
  \begin{eqnarray}\label{3b}
   R_{\mu,\nu}(x,y,z)&=&\int_0^\infty J_{\nu}(xt)J_{\nu}(yt)J_{\mu}(zt)t^{1-\mu} dt\nonumber\\ &=&
    \left\{
      \begin{array}{ll}
        0,  \qquad\qquad \qquad\qquad\qquad\qquad \;\hbox{$z< |x-y|$;} \\
         \frac{(xy)^{\mu-1}\sin ^{\mu-\frac{1}{2}}\theta}{\sqrt{2\pi}z^{\mu}} P^{\frac{1}{2}-\mu}_{\nu-\frac{1}{2}}(\cos \theta),
\; \hbox{ \quad   $|x-y| < z< x+y$;}  \\
     \frac{e^{(\mu-\frac{1}{2})\pi i}\sin((\nu-\mu)\pi)(xy)^{\mu-1}\sinh ^{\mu-\frac{1}{2}} \theta }{(\frac{1}{2}\pi^3)^\frac{1}{2}z^\mu}
Q^{\frac{1}{2}-\mu}_{\nu-\frac{1}{2}}(\cosh \theta),   \; \quad \hbox{$ x+y<z$,}
      \end{array}
    \right.
  \end{eqnarray}
provided  $Re(\mu)>-\frac{1}{2}$, $Re(\nu)>-\frac{1}{2}$,
 and  where here we write
$x^2+y^2-z^2=2xy\cos \theta  $ if $|x-y| < z< x+y$
 and $   z^2-x^2-y^2=2xy\cosh \theta$ if $  x+y < z$ .
The associated  Legendre functions $P^\mu_\nu$ and $Q^\mu_\nu$   are  given in term of hypergeometric function by (see \cite{B}, p.122)
\begin{equation}\label{pppp}
 P^\mu_\nu(x)=\frac{1}{\Gamma(1-\mu)}\left(\frac{1+x}{1-x}\right)^{\frac{\mu}{2}}\; _2F_1\left(\nu+1,-\nu,1-\mu,\;\frac{1-x}{2}\right),\quad
-1<x\leq  1
\end{equation}
and
 \begin{equation}\label{QQQ}
 Q^\mu_\nu(x)= e^{\mu \pi i} \frac{\sqrt{\pi}\Gamma(\mu+\nu+1)(x^2-1)^{\frac{\mu}{2}}}{2^{\nu+1}x^{\mu+\nu+1}\Gamma(\nu+\frac{3}{2})}
 \;_2F_1\left(\frac{\mu+\nu}{2}+1,\frac{\mu+\nu+1}{2},\nu+\frac{3}{2},\;\frac{1}{x^2}\right), 1<x.
 \end{equation}
  It will be observed that  if  $\nu-\mu=n$  is a nonnegative  integer then
\begin{eqnarray*}
     R_{\mu,\nu}(x,y,z)&=&
    \left\{
      \begin{array}{ll}
\frac{2^{\frac{1}{2}-\mu}\Gamma(2\mu) n!}{\Gamma(\nu+\mu)\Gamma(\mu+\frac{1}{2})}\frac{(xy)^{\mu-1}\sin ^{2\mu-1}\phi}{\sqrt{2\pi}z^{\mu}}  C_n^\mu(\cos \theta),
\; \hbox{ $|x-y| < z< x+y$;} \\
        0,  \qquad\qquad \qquad\qquad\qquad\qquad\qquad\;\hbox{$z< |x-y|$ or $z> x+y$.} \\
          \end{array}
    \right.
\end{eqnarray*}
where $C_n^\nu$ is the Gegenbauer polynomial.
\par We shall now discuss   integral representations
which are to be associated with the  Hankel transform. It is a well-known fact from  the theory of  Hankel transform  (see  \cite{Tit}, Ch.8)  that if $f$ is an
integrable function on  $(0,+\infty)$  and of bounded variation in a neighborhood of $t>0$, then the following holds
$$ \int_0^{+\infty}\left\{\int_0^{+\infty}f(r)J_\alpha(rz)\sqrt{rz}\;dr\right\}J_\alpha(tz)\sqrt{tz}\;dz=\frac{f(t+0)+f(t-0)}{2},$$
where $\alpha> -\frac{1}{2}$. If  we take $\alpha=\mu$ and
$$ f(r)= J_{\nu}(xr)J_{\nu}(yr) r^{\frac{1}{2} -\mu}$$
with  $\nu> -\frac{1}{2}$ and  $ \frac{1}{2} <\mu< 2\nu+\frac{3}{2}$ ( which assert the integrability of $f$ ) then  we have
\begin{equation*}
  J_{\nu}(xt)J_{\nu}(yt) t^{ -\mu}= \int_{0}^{\infty}  R_{\mu,\nu}(x,y,z)J_{\mu}(zt)z dz.
\end{equation*}
The formula can be extended  to  $\mu> -\frac{1}{2}$ and $\nu> -\frac{1}{2}$ by the principle of analytic continuation.
Hence in view of (\ref{Jj}) it follows that
\begin{equation}\label{1}
 (xy)^\nu  t^{ 2(\nu-\mu)} \mathcal{J}_{\nu}(xt)\mathcal{J}_{\nu}(yt)= \frac{2^{2\nu-\mu}\Gamma^2(\nu+1)}{ \Gamma(\mu+1)} \int_{0}^{\infty}  R_{\mu,\nu}(x,y,z)\mathcal{J}_{\mu}(zt) z^{\mu+1} dz.
\end{equation}
Taking   $\alpha=\nu$ and
$$ f(r)= J_{\nu}(xr)J_{\mu}(yr) r^{\frac{1}{2} -\mu},$$
a similar argument proves  that
\begin{equation*}
  J_{\nu}(xt)J_{\mu}(yt) t^{ -\mu}= \int_{0}^{\infty}  R_{\mu,\nu}(x,z,y)J_{\nu}(zt)z dz.
\end{equation*}
with   $\nu> -\frac{1}{2}$ and  $\mu> -\frac{1}{2}$. From which we have
\begin{equation}\label{2}
 x^\nu y^\mu   \mathcal{J}_{\nu}(xt)\mathcal{J}_{\mu}(yt)=  2^{\mu}\Gamma(\mu+1)  \int_{0}^{\infty}  R_{\mu,\nu}(x,z,y)\mathcal{J}_{\nu}(zt) z^{\nu+1} dz.
\end{equation}
\par
Let us now   consider the product   $B_{k,a}(\lambda,x)B_{k,a}(\lambda,y)$  which  in virtue of (\ref{BB}) is equal to
\begin{eqnarray}\label{PP}
 && \mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda x|^{a/2}\right) \mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda y|^{a/2}\right)\nonumber
\\&&\qquad\qquad\qquad\qquad\qquad+m_{k,a}^2  \lambda^2xy\mathcal{J}_{\frac{2k+1}{a}}\left(\frac{2}{a}|\lambda x|^{a/2}\right)\mathcal{J}_{\frac{2k+1}{a}}\left(\frac{2}{a}|\lambda y|^{a/2}\right)\nonumber
\\&&+m_{k,a} \lambda x  \mathcal{J}_{\frac{2k+1}{a}}\left(\frac{2}{a}|\lambda x|^{a/2}\right)
 \mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda y|^{a/2}\right)\nonumber
 \\&&\qquad\qquad\qquad \qquad\qquad+m_{k,a}  \lambda y
 \mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda x|^{a/2}\right)\mathcal{J}_{\frac{2k+1}{a}}\left(\frac{2}{a}|\lambda y|^{a/2}\right).
\end{eqnarray}
If we make use (\ref{1}) with $\mu=\nu=\frac{2k-1}{a} $ and $t=\frac{2}{a}|\lambda|^{\frac{a}{2}}$  the first  product of two Bessel functions in (\ref{PP})   may    be  written  as (for $x\neq 0$, $y\neq 0$)
\begin{eqnarray*}
 &&\mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda x|^{a/2}\right) \mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda y|^{a/2}\right)
\\&& = \frac{2^{ \frac{2k-1}{a} }\Gamma(\frac{2k-1}{a}+1)}{|xy|^{k-\frac{1}{2}}}\int_0^\infty R_{\frac{2k-1}{a},\frac{2k-1}{a}}(|x|^{\frac{a}{2}},|y|^{\frac{a}{2}},z)\mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda| ^{a/2}z\right)z^{\frac{2k-1}{a}+1}\; dz
\\ &&= a2^{ \frac{2k-1}{a}-1 }\Gamma\left(\frac{2k-1}{a}+1\right)\int_0^\infty \frac{R_{\frac{2k-1}{a},\frac{2k-1}{a}}(|x|^{\frac{a}{2}},|y|^{\frac{a}{2}},z^{\frac{a}{2}})}{(|xy|z)^{k-\frac{1}{2}}}\mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda| ^{a/2}z^{\frac{a}{2}}\right)z^{ 2k+a-2}\; dz
\\&&= a2^{ \frac{2k-1}{a}-2}\Gamma\left(\frac{2k-1}{a}+1\right)\int_{-\infty}^\infty \frac{R_{\frac{2k-1}{a},\frac{2k-1}{a}}(|x|^{\frac{a}{2}},|y|^{\frac{a}{2}},|z|^{\frac{a}{2}})}{|xyz|^{k-\frac{1}{2}}} B_{k,a}(\lambda,z) |z|^{ 2k+a-2}\; dz.
\end{eqnarray*}
Using  (\ref{1})   with  $\nu=\frac{2k+1}{a} $ and $\mu= \frac{2k-1}{a} $  the second product in (\ref{PP}) can also be written  as
\begin{eqnarray*}
 &&  m_{k,a}^2 \lambda^{ 2}xy \mathcal{J}_{\frac{2k+1}{a}}\left(\frac{2}{a}|\lambda x|^{a/2}\right) \mathcal{J}_{\frac{2k+1}{a}}\left(\frac{2}{a}|\lambda y|^{a/2}\right)=m_{k,a}^2\frac{2^{ \frac{2k-1}{a}} a^{\frac{4}{a}}\Gamma^2( \frac{2k+1}{a}+1)}{\Gamma( \frac{2k-1}{a}+1)}\\&&\qquad\qquad\quad\times   \int_0^{+\infty} sgn(xy)\frac{R_{\frac{2k-1}{a},\frac{2k+1}{a}}(|x|^{\frac{a}{2}},|y|^{\frac{a}{2}},z^{\frac{a}{2}})}{|xyz|^{k-\frac{1}{2}}}\mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda| ^{a/2}z^{\frac{a}{2}}\right)z^{ 2k+a-2}\; dz
\\&&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad= e^{\frac{-2i\pi}{a}}a 2^{ \frac{2k-1}{a}-2}\Gamma\left(\frac{2k-1}{a}+1\right)\\&&\qquad\qquad\qquad\qquad\times\int_{-\infty}^{+\infty} sgn(xy)\frac{R_{\frac{2k+1}{a},\frac{2k+1}{a}}(|x|^{\frac{a}{2}},|y|^{\frac{a}{2}},|z|^{\frac{a}{2}})}{|xyz|^{k-\frac{1}{2}}} B_{k,a}(\lambda,z) |z|^{ 2k+a-2}\; dz.
\end{eqnarray*}
 Applying now in the same manner (\ref2) with $v=\frac{2k+1}{a}$ and $\mu= \frac{2k-1}{a}$ we obtain  that
\begin{eqnarray*}
 && m_{k,a} \lambda x \mathcal{J}_{\frac{2k+1}{a}}\left(\frac{2}{a}|\lambda x|^{a/2}\right) \mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda y|^{a/2}\right)=
  a 2^{ \frac{2k-1}{a}-1}\Gamma\left(\frac{2k-1}{a}+1\right) m_{k,a} \\&&\qquad\qquad\times \int_0^{+\infty} sgn(x)\frac{R_{\frac{2k-1}{a},\frac{2k+1}{a}}(|x|^{\frac{a}{2}},|z|^{\frac{a}{2}},|y|^{\frac{a}{2}})}{(|xy|z)^{k-\frac{1}{2}}}\lambda z\mathcal{J}_{\frac{2k+1}{a}}\left(\frac{2}{a}|\lambda| ^{a/2}z\right)z^{ 2k+a-2}\; dz
\\&& \qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\quad=a 2^{ \frac{2k-1}{a}-2}\Gamma\left(\frac{2k-1}{a}+1\right)\\&&
\qquad\qquad\qquad\qquad\quad\times \int_{-\infty}^{+\infty} \frac{sgn(xz)R_{\frac{2k-1}{a},\frac{2k+1}{a}}(|x|^{\frac{a}{2}},|z|^{\frac{a}{2}},|y|^{\frac{a}{2}})}{|xyz|^{k-\frac{1}{2}}} B_{k,a}(\lambda,z)
 |z|^{ 2k+a-2}\; dz
\end{eqnarray*}
 and
\begin{eqnarray*}
 && m_{k,a} |\lambda|y \mathcal{J}_{\frac{2k+1}{a}}\left(\frac{2}{a}|\lambda y|^{a/2}\right) \mathcal{J}_{\frac{2k-1}{a}}\left(\frac{2}{a}|\lambda x|^{a/2}\right)
 = a 2^{ \frac{2k-1}{a}-2}\Gamma\left(\frac{2k-1}{a}+1\right)\\&&\qquad\qquad\qquad\times\int_{-\infty}^{+\infty} sgn(yz)\frac{R_{\frac{2k-1}{a},\frac{2k-1}{a}}(|y|^{\frac{a}{2}},|z|^{\frac{a}{2}},|x|^{\frac{a}{2}})}{|xyz|^{k-\frac{1}{2}}} B_{k,a}(\lambda,z) |z|^{ 2k+a-2}\; dz.
\end{eqnarray*}
We are thus led to the  formula
\begin{equation}\label{pr}
   B_{k,a}(\lambda,x) B_{k,a}(\lambda,y)= \int_{-\infty}^{+\infty} B_{k,a}(\lambda,z)\Delta_{k,a}(x,y,z) |z|^{ 2k+a-2}\; dz
\end{equation}
 where
\begin{eqnarray*}\label{12}
&&  \Delta_{k,a}(x,y,z)= a 2^{ \frac{2k-1}{a}-2}\Gamma\left(\frac{2k-1}{a}+1\right)\\&&\times \Bigg\{ \frac{R_{\frac{2k-1}{a},\frac{2k-1}{a}}(|x|^{\frac{a}{2}},|y|^{\frac{a}{2}},|z|^{\frac{a}{2}})}{|xyz|^{k-\frac{1}{2}}}
+   e^{\frac{-2i\pi}{a}}sgn(xy)\frac{R_{\frac{2k-1}{a},\frac{2k+1}{a}}(|x|^{\frac{a}{2}},|y|^{\frac{a}{2}},|z|^{\frac{a}{2}})}{|xyz|^{k-\frac{1}{2}}}
  \\&& + sgn(xz)\frac{R_{\frac{2k-1}{a},\frac{2k+1}{a}}(|x|^{\frac{a}{2}},|z|^{\frac{a}{2}},|y|^{\frac{a}{2}})}{|xyz|^{k-\frac{1}{2}}}
 +sgn(yz)\frac{R_{\frac{2k-1}{a},\frac{2k+1}{a}}(|y|^{\frac{a}{2}},|z|^{\frac{a}{2}},|x|^{\frac{a}{2}})}{|xyz|^{k-\frac{1}{2}}}\Bigg\}.
\end{eqnarray*}
\begin{lem}
Let  $\mu>-\frac{1}{2}$ and $\nu>-\frac{1}{2}$. As variables $x>0$ and $y>0$ the integral
 $$\int_{0}^{+\infty}\frac{|R_{ \mu,\nu}( x,y, z)|}{(xy)^{ \mu}}\;z^{\mu+1}dz$$
is uniformly bounded.
\end{lem}
\begin{proof}
  The proof is based on  the integrals  of  \cite{WM1}    that  appeared in  (16) of 18.1 and  (23) and of 18.2, to get the following
\begin{eqnarray}\label{P}
 \int_{-1}^{1}(1-t^2)^{\frac{\mu}{2}-\frac{1}{4}}P_{\nu-\frac{1}{2}}^{\frac{1}{2}-\mu}(t)\;dt= \frac{\pi2^{\frac{1}{2}-\mu}\Gamma(\mu+ \frac{1}{2}) }{\left(\Gamma(\frac{\mu+\nu+1}{2})\right)^2 \Gamma(\frac{\mu-\nu+2}{2}) \Gamma(\frac{\mu-\nu+1}{2}) },
\end{eqnarray}
and
\begin{eqnarray}\label{Q}
 \int_{1}^+{\infty}(t^2-1)^{\frac{\mu}{2}-\frac{1}{4}}Q_{\nu-\frac{1}{2}}^{\frac{1}{2}-\mu}(t)\;dt=\sqrt{2} e^{i(\frac{1}{2}-\mu\pi)}
 \frac{\Gamma(\frac{1+\nu-\mu}{2})\Gamma(\frac{\nu-\mu}{2}+\frac{1}{4})\Gamma(\mu+\frac{3}{4})\Gamma(\frac{3}{4})}{\Gamma(\nu+\mu)\Gamma(\nu+\mu+1)}.
\end{eqnarray}
From (\ref{3b}) we have
 \begin{eqnarray*}
\int_{x+y}^{+\infty}\frac{|R_{ \mu,\nu}( x,y, z)|}{(xy)^{ \mu}}z^{\mu+1}dz=  \frac{ |\sin((\nu-\mu)\pi)|}{(\frac{1}{2}\pi^3)^\frac{1}{2}} \int_{x+y}^{+\infty}\frac{ \sinh ^{\mu-\frac{1}{2}} \theta }{xy }
Q^{\frac{1}{2}-\mu}_{\nu-\frac{1}{2}}(\cosh \theta)z dz.
\end{eqnarray*}
Putting the change of variable
$$t=\cosh \theta=\frac{z^2-x^2-y^2}{2xy}, $$
it follows that
\begin{eqnarray}\label{QQ}
\int_{x+y}^{+\infty}\frac{|R_{ \mu,\nu}( x,y, z)|}{(xy)^{ \mu}}z^{\mu+1}dz=  \frac{ |\sin((\nu-\mu)\pi)|}{(\frac{1}{2}\pi^3)^\frac{1}{2}} \int_{1}^{+\infty}
 (t^2-1)^{\frac{\mu}{2}-\frac{1}{4}}Q_{\nu-\frac{1}{2}}^{\frac{1}{2}-\mu}(t)\;dt.
\end{eqnarray}
 Similarly
  \begin{eqnarray*}
\int_{|x-y|}^{x+y}\frac{|R_{ \mu,\nu}( x,y, z)|}{(xy)^{ \mu}}z^{\mu+1}dz= \frac{1}{\sqrt{2\pi}} \int_{-1}^{1} (1-t^2)^{\frac{\mu}{2}-\frac{1}{4}}
 | P^{\frac{1}{2}-\mu}_{\nu-\frac{1}{2}}( t )|\;dt.
\end{eqnarray*}
In view of (\ref{pppp}) we see that    $P^{\frac{1}{2}-\mu}_{\nu-\frac{1}{2}}( t )\geq 0$ when $-\frac{1}{2}<\nu\leq \frac{1}{2}$ . Thus using (\ref{P}) together  with the contiguous relation (see 4.3.3 of \cite{WM}),
$$ P_{\nu+1}^\mu(t)= tP_{\nu}^\mu(t)-(\mu+\nu)(1-t^2)^{\frac{1}{2}}P_{\nu}^{\mu-1}(t)$$
 one can see that
$$\int_{|x-y|}^{x+y}\frac{|R_{ \mu,\nu}( x,y, z)|}{(xy)^{ \mu}}z^{\mu+1}dz$$
is uniformly bounded. Then combine this with (\ref{QQ}) and (\ref{Q}) to achive  the proof of the lemma.
\end{proof}
\begin{lem}
For $\mu>-\frac{1}{2}$ and $\nu>-\frac{1}{2}$ the integral
 $$\int_{0}^{+\infty}\frac{|R_{ \mu,\nu}( x,z, y)|}{(xy)^{ \mu}}\;z^{\mu+1}dz$$
is uniformly bounded  with respect to $x>0$ and $y>0$.
\end{lem}
\begin{proof} Let us denote by
$$I_1(x,y)=\int_{|x-y|}^{x+y}\frac{|R_{ \mu,\nu}( x,z, y)|}{(xy)^{ \mu}}\;z^{\mu+1}dz
\quad\text{and}\quad I_2(x,y)=\int_{x+y}^{\infty}\frac{|R_{ \mu,\nu}( x,z, y)|}{(xy)^{ \mu}}\;z^{\mu+1}dz.$$
We are therefore led to prove that $I_1(x,y)$ and $I_2(x,y)$ are bounded.
It is convenient to divide the proof into two cases  $x\geq y$  and   $x < y$. We use the letter $C$   to denote positive constant whose value  can change at each occurrence.
 \par Let us begin with the case $x\geq y$ where  we have  $I_2(x,y)=0$.  To establish the boundedness of $I_1$   we use the following identity
  \begin{eqnarray}
    \Gamma(1-\mu) P_{\nu}^\mu(t)=2^\mu(1-t^2)^{-\frac{\mu}{2}}\;_2F_1 \left( \frac{1+\nu-\mu}{2}, \frac{-\mu-\nu}{2}, 1-\mu,1-t^2 \right)
  \end{eqnarray}
   which follows   from well known properties of the hypergeometric function  $_2F_1$ (see also \cite{WM}, p.167).
In addition  the function $\;_2F_1 \left( \frac{1+\nu-\mu}{2}, \frac{-\mu-\nu}{2}, 1-\mu,1-t^2 \right)$ is bounded
when $0<t<1$. It is then clear that
\begin{equation}\label{ppp}
  |P_{\nu}^\mu(t)|\leq C \;(1-t^2)^{-\frac{\mu}{2}},\qquad 0\leq t\leq 1.
\end{equation}
Now  using (\ref{ppp}), we get
when $|x-z|\leq y \leq x+z$ ( which is also equivalent to  $x-y\leq z \leq x+y$),
$$\frac{|R_{ \mu,\nu}( x,z, y)|}{(xy)^\mu}\leq C\; \frac{z^{\mu-1}}{  x y^{2\mu}} \left\{1- \left(\frac{x^2+z^2-y^2}{2xz}\right)^2\right\}^{\mu-\frac{1}{2}}.$$
For convenience, we write
 $$ 1- \left(\frac{x^2+z^2-y^2}{2xz}\right)^2= \frac{((x+y)^2-z^2)(z^2-(x-y)^2)}{4(xz)^2}.$$
  Hence,
  $$\frac{|R_{ \mu,\nu}( x,z, y)|}{(xy)^\mu}\leq C \; \frac{\Big\{((x+y)^2-z^2)(z^2-(x-y)^2\Big\}^{\mu-\frac{1}{2}}}{(xyz)^{2\mu}}\;z^\mu=CW(x,y,z)z^{\mu}.$$
  Now observe that
  $$ \int_{x-y}^{x+y}W(x,y,z)z^{2\mu+1}\;dz=\frac{2^{2\mu-1}\sqrt{\pi}\Gamma(\mu+\frac{1}{2})}{\Gamma(\mu+1)}$$
and therefore we conclude that  $I_1(x,y)$ is   bounded.  Consider now    $y\geq x$.   We shall use  the following estimates that follows from   (\ref{pppp})  and  15.4(ii) of  \cite{Nist},
\begin{eqnarray}
|P_\nu^\mu(t)|&\leq &C (1-t^2)^{-\frac{\mu}{2} }, \qquad\text{if }\qquad     \mu >0,\label{11}\\
 |P_\nu^\mu(t)|&\leq& C (1-t^2)^{\frac{\mu}{2}},\qquad\text{if }\qquad   \mu <0,\label{22}\\
|P_\nu^\mu(t)|&\leq& C |\ln(e(1+t))| , \qquad\text{if } \qquad     \mu =0,\label{33}
\end{eqnarray}
where $-1<t<1$. Noting first that in view of (\ref{11}) and  (\ref{ppp}) one can  conclude the boundedness of $I_1$ for  $ \mu <\frac{1}{2}$   in a similar manner as before.
When $ \mu > \frac{1}{2}$  and  from (\ref{22}) we have  for $y-x<z<x+y$,
$$\frac{|R_{ \mu,\nu}( x,z, y)|}{(xy)^\mu}\leq C \;\frac{z^{\mu-1}}{xy^{2\mu}}$$
 and  thus,
 $$I_1(x,y)=\int_{y-x}^{x+y}\frac{|R_{ \mu,\nu}( x,z, y)|}{(xy)^{ \mu}}\;z^{\mu+1}dz\leq C\; \frac{(x+y)^{2\mu+1}-(y-x)^{2\mu+1}}{xy^{2\mu}}$$
 $$\leq C\; \frac{(x/y+1)^{2\mu+1}-(1-x/y)^{2\mu+1}}{ x/y}\leq C.$$
Since the function $\Big(t+1)^{2\mu+1}-(1-t)^{2\mu+1}\Big)t^{-1}$ is bounded on $(0,1)$. In the case $\mu=\frac{1}{2}$, the estimation of (\ref{33}) gives $$|I_1(x,y)|\leq \frac{C}{xy}\int_{y-x}^{x+y}\left(1+\ln\left(1+\frac{x^2+z^2-y^2}{2xz}\right)\right)zdz$$
 Using the Change of  variable
$$t=\frac{x^2+z^2-y^2}{2xz},$$
 one can see that
$$\frac{1}{xy}\int_{y-x}^{x+y}\ln\left(1+\frac{x^2+z^2-y^2}{2xz}\right)zdz\leq2 \int_{-1}^1 \frac{\ln(1+t)}{|t|}\;dt.$$
As a consequence $I_1$ is bounded.
We  come now to the boundedness of $I_2$.  According with (\ref{QQQ}) and  15.4(ii) of  \cite{Nist} we  get
\begin{eqnarray}
 |Q_\nu^\mu(t)|&\leq& C \frac{(t^2-1)^{-\frac{\mu}{2} }}{t^{\nu-\mu+1}},,\qquad\text{if }\qquad   \mu >0,\label{55}\\
|Q_\nu^\mu(t)|&\leq &C \frac{(t^2-1)^{\frac{\mu}{2} }}{t^{\nu+\mu+1}}, \qquad\text{if }\qquad     \mu <0,\label{44}\\
|Q_\nu^\mu(t)|&\leq& C \frac{(t^2-1)^{\frac{\mu}{2} }}{t^{\mu+\nu+1}}|\ln(1-t^{-2})| , \qquad\text{if } \qquad     \mu =0. \label{66}
\end{eqnarray}
If $\mu  > \frac{1}{2}$ then under  consideration  (\ref{44})   with  (\ref{3b}) we have
$$\frac{|R_{ \mu,\nu}( x,z, y)|}{(xy)^\mu}\leq C x^{\nu-\mu}y^{-2\mu}(y^2-x^2-z^2)^{\mu-\nu-1}z^{\nu}$$
 and
\begin{eqnarray*}
| I_2(x,y)|&\leq C& x^{\nu-\mu}y^{-2\mu} \int_{0}^{y-x} \frac{z^{\mu+\nu+1}}{(y^2-x^2-z^2)^{\nu-\mu+1}}\; dz
 \\&\leq&Cx^{\nu-\mu}y^{-2\mu}(y^2-x^2)^{\frac{3\mu-\nu}{2}} \int_{0}^{ \sqrt{\frac{y-x}{y+x}}} \frac{z^{\mu+\nu+1}}{(1-z^2)^{\nu-\mu+1}}\; dz
 \\&\leq& C \Psi(x/y),
\end{eqnarray*}
where
$$\Psi(t)=t^{\nu-\mu}(1-t^2)^{\frac{3\mu-\nu}{2}} \int_{0}^{ \sqrt{\frac{1-t}{1+t}}} \frac{z^{\mu+\nu+1}}{(1-z^2)^{\nu-\mu+1}}\; dz.$$
It  not hard  to verify that $\Psi$ is bounded on $(0,1)$, which implies that $I_2$ is bounded.\\
If $\mu<\frac{1}{2}$ then
$$ |I_2(x,y)|\leq \frac{C}{xy^{2\mu}}\int_0^{y-x}\frac{\left\{\left(\frac{y^2-x^2-z^2}{2xz}\right)^2-1\right\}^{\mu-\frac{1}{2}}}{\left(\frac{y^2-x^2-z^2}{2xz}\right)^{\nu+\mu}}
z^{2\mu}dz$$
letting  the  change of variable
$$t=\frac{y^2-x^2-z^2}{2xz},$$
it becomes
$$ |I_2(x,y)|\leq C\;y^{-2\mu}\int_1^{+\infty } \frac{(t^2-1)^{\mu-\frac{1}{2}}}{t^{\nu+\mu}}\;\frac{(\sqrt{x^2t^2+y^2-x^2}-xt)^{2\mu+1}}{\sqrt{x^2t^2+y^2-x^2}}\; dt .$$
As $y>x$
  $$  \frac{(\sqrt{x^2t^2+y^2-x^2}-xt)^{2\mu+1}}{\sqrt{x^2t^2+y^2-x^2}}\leq \left(\frac{y^2-x^2}{y} \right)^{2\mu+1}\leq y^{2\mu},$$
it follows that
$$ |I_2(x,y)|\leq C\; \int_1^{+\infty} \frac{(t^2-1)^{\mu-\frac{1}{2}}}{t^{\nu+\mu}} \; dt .$$
Similarly,  when  $\mu=\frac{1}{2}$ where  it follows from (\ref{66}) that
$$ |I_2(x,y)|\leq C\; \int_1^{+\infty} \frac{ \ln(1-t^{-2})}{t^{\nu+1/2}} \; dt .$$
  Consequently, the boundedness  of $I_2$ follows. This completes the proof of the lemma.
  \par Now our main result can be stated as follows.
\begin{thm}
In one dimentional case the kernel
$ B_{k,a}$ satisfies  the product formula
$$ B_{k,a}(\lambda,x) B_{k,a}(\lambda,y)= \int_{-\infty}^{+\infty} B_{k,a}(\lambda,z)d\gamma_{x,y}^{k,a}(z)$$
where
$$d\gamma_{x,y}^{k,a}(z)=\left\{
  \begin{array}{ll}
     \Delta_{k,a}(x,y,z) |z|^{ 2k+a-2}dz, & \hbox{if $xy\neq 0$;} \\
   \delta_x (z) , & \hbox{if $y=0$  ;} \\
     \delta_y(z) & \hbox{if $x=0$.}
  \end{array}
\right.$$
Further the integral $$ \int_{-\infty}^{+\infty}|d\gamma_{x,y}^{k,a}(z) |$$
   is uniformly bounded.
\end{thm}
  Noting here  that the measure $\delta_{x,y}^{k,a}$  has compact  support  if and only if $a=\frac{2}{n}$, $n\in \mathbb{N}$.
  \end{proof}