I've only recently read up on this stuff, so take everything here with a grain of salt. What we can do is try to mimic the idea that the Laplacian is the divergence of the gradient. So, we first define these two operations.
Defining "the" Gradient.
First, we address the gradient. Given a $(p,q)$ tensor field $T$, the covariant derivative $\nabla T$ will be a $(p, q+1)$ tensor field. Using the metric $g$'s musical isomorphism, we can construct a $(p+1,q)$ tensor field. For the lack of a more creative name, I'll call it $\text{grad}(T)$. So,
\begin{align}
[\text{grad}(T)](\omega_1, \dots \omega_p, \omega, X_1, \dots, X_q) &:= \left(\nabla_{g^{\sharp}(\omega)}T\right)(\omega_1, \dots, \omega_p, X_1, \dots, X_q)
\end{align}
where the $\omega$'s are covector fields, and $X$'s are vector fields. Or in components,
\begin{align}
[\text{grad}(T)]^{i_1 \dots i_p, a}_{j_1 \dots j_q} &= \left(\nabla_{g^{as}\frac{\partial}{\partial x^s}}T \right)^{i_1 \dots i_p}_{j_1 \dots j_q} = g^{as}\left(\nabla_sT \right)^{i_1 \dots i_p}_{j_1 \dots j_q}.
\end{align}
In the case of a smooth function $f \in C^{\infty}(M)$ this gives us that $\text{grad}(f) = g^{\sharp}(df)$, so it is indeed the vector field "associated" to $df$ via the metric.
Defining "the" Divergence.
Next, we consider the divergence. Let $T$ be any $(p,q)$ tensor field with $p \geq 1$. Then, $\nabla T$ will be a $(p, q+1)$ tensor field, and we can take a trace with respect to the last covariant and contravariant slots to end up with a $(p-1, q)$ tensor field, which we call $\text{div}(T)$:
\begin{align}
[\text{div}(T)](\omega_1, \dots \omega_{p-1}, X_1, \dots, X_q) &:= \text{trace} \bigg((\omega,X) \mapsto (\nabla T) \left(\omega_1, \dots, \omega_{p-1}, \omega, X_1, \dots, X_q, X \right) \bigg)
\end{align}
where the trace of a $(1,1)$ tensor (field) is of course well-defined and chart/basis/coordinate independent. So, in components, it is given by
\begin{align}
[\text{div}(T)]^{i_1 \dots i_{p-1}}_{j_1 \dots j_q} &= (\nabla_sT)^{i_1 \dots i_{p-1},s}_{j_1 \dots j_q}
\end{align}
Defining "the" Laplacian
So, we can now define a Laplacian of any $(p,q)$ tensor field by:
\begin{align}
\text{Lap}(T) &:= \text{div}(\text{grad}(T)).
\end{align}
This makes sense because $\text{grad}(T)$ will be a $(p +1, q)$ tensor field, so we can take its divergence to get a $(p,q)$ tensor field again.
For smooth functions, this definition of the Laplacian yields, in coordinates,
\begin{align}
\text{Lap}(f) &= \dfrac{1}{\sqrt{|g|}} \dfrac{\partial}{\partial x^a} \left( \sqrt{|g|}\, g^{as} \dfrac{\partial f}{\partial x^s}\right),
\end{align}
and assuming I didn't make any sign mistakes in my calculation, I'm pretty sure it also coincides with the definition you gave: $\star d (\star df)$.
For the definition of divergence, I used the one presented in Sachs and Wu's General Relativity book. Here are some observations I made: the definition of gradient required us to "raise" one of the indices, and I simply chose the last one. Similarly, the definition of divergence requires us to choose a certain pair of covariant and contravariant indices to take the trace over; once again I chose (as in the book) the last two. If the tensor $T$ is symmetric, then it of course doesn't matter which ones we choose. If it is not symmetric, then I haven't really carefully considered how it affects the definitions.
Finally, I'm no expert in this, so I'm not sure how exactly this relates to other definitions you might find.
