On the recurrence $a_{k+1}=\frac{a_k}{p_{k+1}}+\frac{1}{a_k}$, with $a_1=1$ and being $p_k$ the $kth$ prime number We denote for integers $n\geq 1$ the $nth$ prime number as $p_n$. And for integers $k\geq 1$, I consider this recurrence relation
$$a_{k+1}=\frac{a_k}{p_{k+1}}+\frac{1}{a_k},\tag{1}$$ with $a_1$ defined to be equal to $1$.
I've calculated some (few) terms of this sequence $a_k$. See here the first  examples.
Examples of computations of some terms of our sequence $ \left\{ a_n\right\}_{n=1} ^\infty$: 
1) Since $a_1=1$ then $a_2=\frac{a_1}{p_2}+(a_1)^{-1}=\frac{1}{3}+\frac{1}{1}=\frac{4}{3}$. 
2) Since the third prime number is $p_3=5$ one has  $$a_3=\frac{a_2}{p_3}+\frac{1}{a_2}=\frac{4/3}{5}+\frac{3}{4}=\frac{61}{60}.$$ 
3) Similarly $a_4=\frac{61}{7\cdot 60}+(60/61)^{-1}\approx 1.1288$. $\square$
Thus our sequence starts as $a_1=1,a_2\approx 1.3333, a_3\approx 1.0117, a_4\approx 1.1288$ and similarly we can calculate $a_5\approx 0.9664,a_6\approx 1.1091$ or $a_7\approx 0.9669$.
To create this problem I was inspired in a recurrence that $\sqrt{2}$ solves as you can see from this WIkipedia, and now I am curious about how to check that the sequence in $(1)$ is convergent.

Question. Please, can you prove that $\left\{ a_n\right\}_{n=1} ^\infty$ defined from $(1)$ is convergent? Many thanks.

I know a a main tool in the theory of prime numbers is the Prime Number Theorem: $p_n\sim n\log n$ as $n\to\infty$.
 A: This partially answers the question by showing that both $(a_{2n})$ and $(a_{2n+1})$ converge. The proof may not be an easy path, and particularly the first section can be quite boring. You may focus on only color-boxed statements and skip the rest of the preliminary.

1. Preliminary

Definition. For $p \geq 2$ we define $f_p : (0, \infty) \to (0, \infty)$ by
$$f_p(x) = \frac{x}{p} + \frac{1}{x}$$
Also we set $g_{p,q} = f_q \circ f_p$ for $q \geq p \geq 2$.

We would like to investigate these functions. Then it is immediate that the equation $f_p'(x) = 0$ has a unique zero $x = \sqrt{p}$ on $(0, \infty)$. Using this, we can solve $g_{p, q}'(x) = f_q'(f_p(x)) f_p'(x) = 0$.


*

*One zero comes from $f_p'(x) = 0$, yielding $x = \sqrt{p}$.

*There are two other zeros that come from $f_q'(f_p(x)) = 0$ or equivalently $f_p(x) = \sqrt{q}$. Solving this, we obtain two zeros
$$ \bbox[#fff9e3,border:1px solid #ffeb8e,12px]{ \alpha_{p,q} := \frac{p\sqrt{q} - \sqrt{ p (pq - 4)} }{2} }
\quad \text{and} \quad
\beta_{p,q} := \frac{p\sqrt{q} + \sqrt{ p (pq - 4)} }{2} \tag{1} $$
As for the smaller zero, we have $\alpha_{p,q} = \frac{2}{\sqrt{q} + \sqrt{q - (4/p)}} < \frac{2}{\sqrt{q}} $ and hence it tends zero as $q \to \infty$. Similarly, we have $\beta_{p,q} \geq p\sqrt{q}/2 \geq \sqrt{p}$, where the last inequality follows from $q \geq p \geq 2$.
Combining these observations, we obtain

Lemma 1. $g_{p, q}$ is strictly increasing on $I_{p,q} := [\alpha_{p,q}, \sqrt{p}]$.

Next we investigate the fixed point of $g_{p,q}$. By a brutal-force computation, we find that $g_{p,q}$ has the unique fixed point
$$ \bbox[#fff9e3,border:1px solid #ffeb8e,12px]{ x_{p,q} := \sqrt{\frac{p}{\sqrt{pq} - 1}} } \tag{2}$$
on $(0, \infty)$. Moreover, since $g_{p,q}(\alpha_{p,q}) = \frac{2}{\sqrt{q}} > \alpha_{p,q}$ and $g_{p,q}(\sqrt{p}) = f_q(2/\sqrt{p}) \leq \sqrt{p}$, IVT tells that $x_{p,q} \in I_{p,q}$. So it follows that

Lemma 2. We have
(i) If $x \in [\alpha_{p,q}, x_{p,q})$, then $x < g_{p,q}(x) < x_{p,q}$. 
(ii) If $x \in (x_{p,q}, \sqrt{p}]$ then $x > g_{p,q}(x) > x_{p,q}$.

This situation is summarized in the following graph:
$\hspace{7em}$ 
We also need a bit technical observation.

Lemma 3. Let $q' \geq p' \geq 2$ be such that $q' \geq q$. Then $g_{p,q}(I_{p,q}) \subseteq I_{p',q'}$.

In order to prove this, it is enough to notice that $ \alpha_{p',q'}
\leq \frac{2}{\sqrt{q'}}
\leq \frac{2}{\sqrt{q}}
= g_{p,q}(\alpha_{p,q}) $
and that
$g_{p,q}(\sqrt{p}) = 2/\sqrt{p} \leq \sqrt{2} \leq \sqrt{p'}$.
Finally we need an input from number theory. Let $p_n$ be the $n$-the smallest prime number. Then

Lemma 4. $\lim_{n\to\infty} p_{n+1} / p_n = 1$.

Proof. It is equivalent to saying that the prime gap satisfies $(p_{n+1} - p_n)/p_n \to 0$. See this for the reference.

2. Main proof
Let $x_n = x_{p_n,p_{n+1}}$, where the right-hand side is the symbol defined as $\text{(2)}$. In a similar fasion, we write $g_n = g_{p_n,p_{n+1}}$ and $I_n = I_{p_n,p_{n+1}}$. Then by Lemma 3, we know that $g_n(I_n) \subseteq I_{n+2}$. Also we have $a_2 = \frac{4}{3} \in [2/\sqrt{5}, \sqrt{3}] = I_2$ and $a_3 = \frac{61}{60} \in [2/\sqrt{7}, \sqrt{5}] = I_3$. So by the induction together with the recurrence relation $a_{n+2} = g_n(a_n)$ tells that $a_n \in I_n$ for all $n \geq 2$.
Now by Lemma 4, we know that $x_n \to 1$ as $n\to\infty$. So if we fix a sufficiently small $\epsilon > 0$, there exists $N = N(\epsilon) \geq 2$ such that $|x_n - 1| < \epsilon$ for all $n \geq N$. So by Lemma 1,


*

*If $a_n \leq 1-\epsilon$, then $a_n < x_n$ and hence $a_n < a_{n+2} < x_n < 1+\epsilon$.

*If $a_n \geq 1+\epsilon$, then $a_n > x_n$ and hence $a_n > a_{n+2} > x_n > 1-\epsilon$.

*Since $a_{n+2}$ always lies between $a_n$ and $x_n$, if $|a_n - 1| < \epsilon$ then $|a_{n+2} - 1| < \epsilon$ as well. 


Now let us fix $r \in \{1,2\}$. Then the above observation tells that we have the following trichotomy for the sequence $(a_{2n+r})_{n=0}^{\infty}$.
Case 1. There exists $\epsilon > 0$ such that $a_{2n+r} \leq 1-\epsilon$ for all large $n$. In this case, the sequence is eventually monotone increasing and hence converges.
Case 2. There exists $\epsilon > 0$ such that $a_{2n+r} \geq 1+\epsilon$ for all large $n$. By a similar reasoning, the sequence converges.
Case 3. For all $\epsilon > 0$ we have $|a_{2n+r} - 1| < \epsilon$ for all large $n$. This immediately translates to the statement that $a_{2n+r}$ converges to $1$.
Combining altogether, we have proved that

Proposition. Both $(a_{2n})$ and $(a_{2n+1})$ converge.

Remark. Calibrating the initial value, this proof should work for any sequence $(p_n)$ such that $p_n \geq 2$ and $p_n \nearrow \infty$ with $p_{n+1}/p_n \to 1$. I strongly suspect that we may find some sequence $(p_n)$ such that all these conditions are met but the limit of $(a_{2n})$ and $(a_{2n+1})$ do not coincide. (Geometric sequence might be such a candidate.) So we indeed need more input to settle down the issue of convergence of $(a_n)$.

3. Numerical computation
Here is a computation of first $10^5$ terms using Mathematica:
$\hspace{4em}$ 
Although it is not easy to read out the convergent bahavior, it seems that the error decays at least as fast as the speed of $1/\log n$. Scaling up the error by log factor indeed provides a better picture:
$\hspace{8em}$ 
A: Download Pari/GP and run this code
 N = 10^5;
 p = primes(N);
 a = vector(N); a[1]=1.0; /* 1.0 means float, with a[1] = 1; instead the sequence of rationals is computed in closed form */
 for( k = 1,  N-1, { a[k+1] = a[k]/p[k+1]+1/a[k];});

 x = vector(N); for( k = 1,  N, { x[k] = k;});
 plothraw(x,a);

A: The answer given in the meantime is very iteresting. However, as you suggested, I write my own results as answer instead of comment, even though I am far from a conclusion. 
I would start from an "easier" situation, with the sequence
$$ a_{k+1} = \frac{a_k}{M} + \frac{1}{a_k},$$
with $M>2$.
I now analyse the function
$$y_M(x) = \frac{x}{M}+\frac{1}{x}$$
and interesect it with $y=x$ to look for the equilibrium point which occurs at
$$x_C= \sqrt\frac{M}{M-1}.$$
The function reaches its minimum at 
$$x_m = \sqrt M.$$
If you take $x=f_M(x_c)=\frac{2\sqrt M}{M}$ and iterate the sequence once more, you get the point
$$
f_M\left(\frac{2\sqrt M}{M}\right)=\sqrt M\left(\frac{2}{M^2}+\frac{1}{2}\right) =x_l
$$
that, for $M>2$ is less then $\sqrt M$. As a consequence,  once a value
$$
\frac{2\sqrt{M}}{M} \leq a_k \leq \sqrt M\left(\frac{2}{M^2}+\frac{1}{2}\right)\tag{1}\label{eq:one}
$$
is reached, the sequence values remain "trapped" in the range
$$ \frac{2\sqrt{M}}{M} <a_{k+1}< \sqrt M$$
and converge to $x_C$. 
So:


*

*When the starting point $a_1$ is in the range 
$$f^{-1}_M(x_l)<a_k<x_l,$$
where the left extremity is taken to be the smallest between the two solutions, then
the sequence will oscillate with decreasing amplitude and converge to $x_C$. 

*If $a_1 > x_l$  then the sequence will first decrease monotonically until the range \eqref{eq:one} is reached; after that the it will behave as in 1.

*Finally, if $a_1 < f^{-1}_M(x_l)$, then from $a_2$ onwards, the sequence will behave as in 2.


As I said, no conclusions here regarding your sequence, in which $M$ is a (growing) function of $k$. The consequences are, I believe:
a. Affecting the ranges at each iteration.
b. As $M$ increases, the sequences approaches the behavior of the sequence $a_k=\frac{1}{a_k}$, which of course does not converge.
That is all for now. 
