A generalization of the product of harmonic numbers to non-integer arguments

This question is somewhat related to one of my previous questions: Fibonorial of a fractional or complex argument.

Recall the definition of harmonic numbers: $$H_n=\sum_{k=1}^n\frac1k=1+\frac12+\,...\,+\frac1n\tag1$$ Obviously, harmonic numbers satisfy the following functional equation: $$H_n-H_{n-1}=\frac1n\tag2$$ The definition $(1)$ is valid only for $n\in\mathbb N$, but it can be generalized to all positive indices. There are several equivalent ways to do this: $$H_a=\sum_{k=1}^\infty\left(\frac1k-\frac1{k+a}\right)=\int_0^1\frac{1-x^a}{1-x}\,dx=\frac{\Gamma'(a+1)}{\Gamma(a+1)}+\gamma\tag3$$ This generalized definition gives a real-analytic function (that can be extended to a complex-analytic if needed) and still satisfies the functional equation $(2)$ even for non-integer values of $a$.

Now, consider the product of harmonic numbers: $$P_n=\prod_{k=1}^nH_k=H_1\,H_2\,...H_n=1\times\left(1+\frac12\right)\times\,...\times\left(1+\frac12+\,...+\frac1n\right)\tag4$$ The numerators and denominators of the terms of this sequence appear as A097423 and A097424 in the OEIS. Obviously, the following function equations hold: $$\frac{P_n}{P_{n-1}}=H_n,\quad\quad\frac{P_n}{P_{n-1}}-\frac{P_{n-1}}{P_{n-2}}=\frac1n\tag5$$ I'm looking for a continuous generalization $P_a$ of the discrete sequence $P_n$, which is real-analytic for all $a>0$ and satisfies the functional equations $(5)$.

Could you suggest a way to construct such a function? Is there a series or integral representation for it? Can we generalize it to complex arguments?

Update: It seems we can use the same trick that is used to define $\Gamma$-function using a limit involving factorials of integers: $$P_a=\lim_{n\to\infty}\left[\left(H_n\right)^a\cdot\prod_{k=1}^n\frac{H_k}{H_{a+k}}\right]=\frac1{H_{a+1}}\cdot\prod_{n=1}^\infty\frac{\left(H_{n+1}\right)^{a+1}}{\left(H_n\right)^a\,H_{a+n+1}}\tag6$$

One may use the Euler-Maclaurin formula to find that as $n\to\infty$, we have

$$\sum_{k=1}^n\ln(H_k)\sim c+\int_0^n\ln(H_x)\ dx+\frac{\ln(H_n)}2+\frac{H_n'}{4H_n}+\dots$$

where $c=\lim_{n\to\infty}\sum_{k=1}^n\ln(H_k)-\left(\int_0^n\ln(H_x)\ dx+\frac{\ln(H_n)}2+\frac{H_n'}{4H_n}\right)\approx0.7288$

One can then use $(5)$ to recursively relate back to the original sum in the same manner you did in $(6)$, as this expansion approximates best as $n\to\infty$.