Asymptotics for the Dirichlet convolution $a * \varphi = 2a - \epsilon$ and the roots of $2\zeta(s) = \zeta(s-1)$

Consider the sequence defined by $a_1 = 1$ and the recurrence relation for $n \ge 2$:

$$a_n = \sum_{k=1}^{n-1} a_{\gcd(n,k)}$$

Grouping the terms by their divisors $d = \gcd(n,k)$, the number of integers $k < n$ such that $\gcd(n,k) = d$ is given by $\varphi(n/d)$, where $\varphi$ is Euler's totient function. This allows the recurrence to be rewritten as $a_n = \sum_{d|n, d<n} a_d \varphi\left(\frac{n}{d}\right)$. Adding $a_n \varphi(1)$ to both sides produces the Dirichlet convolution $2a = a * \varphi$ valid for $n \ge 2$. Over all positive integers, we obtain the identity $a * \varphi = 2a - \epsilon$, where $\epsilon$ is the Dirichlet identity.

Let $A(s) = \sum_{n=1}^\infty a_n n^{-s}$ denote the Dirichlet series for $a_n$. Using the standard generating function $\sum_{n=1}^\infty \varphi(n) n^{-s} = \frac{\zeta(s-1)}{\zeta(s)}$, we obtain the relation $A(s) \frac{\zeta(s-1)}{\zeta(s)} = 2A(s) - 1$. Solving for $A(s)$ yields:

$$A(s) = \frac{\zeta(s)}{2\zeta(s) - \zeta(s-1)}$$

Evaluating the sequence at prime indices yields $a_p = p-1$, establishing the lower limit $\liminf_{n \to \infty} \frac{a_n}{n} = 1$. The maximal order is highly composite, peaking at primorials $N_k = \prod_{i=1}^k p_i$. The dominant real pole of $A(s)$ occurs at the root of $2\zeta(s) = \zeta(s-1)$ in the interval $(2,3)$, which evaluates numerically to $\rho \approx 2.4841$.

My Attempt: To establish the average order of the summatory function $S(x) = \sum_{n \le x} a_n$, we apply Perron's formula. By shifting the contour of integration past the vertical line $\Re(s) = \rho$, we extract the main term $S(x) = c x^\rho + O(E(x))$. However, completing the contour integration becomes intractable when bounding the error term $E(x)$, as this requires controlling the complex roots of the denominator $2\zeta(s) - \zeta(s-1)$. Standard methods relying on the zero-free region of the Riemann zeta function do not trivially extend to linear combinations of shifted zeta functions, leaving the exact contour shift and subsequent error estimation unproven.

Is there established literature on the zero-free regions for linear combinations of shifted zeta functions of the form $c_1\zeta(s) + c_2\zeta(s-k)$ that can be applied to rigorously bound the error term in this contour integration?