There are divisibility sequences of the form $\frac{p^k-1}{p-1}$ which have the property that if $a \mid b$, then $f(a) \mid f(b)$, ensuring (among other things) that only prime indices of the sequence can possibly be prime. This form grows as $2^n$ with $p=2$.
You can find significantly slower growth with the Fibonacci sequence, where the $n$th term is $ \frac{\varphi^n-\psi^n}{\sqrt 5} \approx 1.6^n$.
The slowest-growing form I know of is the recurrence relation defined by $S_n = S_{n-1} -2 S_{n-2}$ with $S_0 = 0, S_1=1$. I don't know the closed-form expression offhand, but I think it grows as something closer to $1.4^n$.
So my question is, are there known slower-growing non-trivial divisibility sequences? By trivial, I'm referring to $\mathbb{N}$ or similar. I would guess that you should be able to generate one with an exponent arbitrarily close to $1$, but that's just speculation. Better yet, of course, would be one with sub-exponential growth, but I wouldn't be surprised if that were impossible.
Loosely related question: it seems like every prime factor $a$ appearing in $\frac{p^k-1}{p-1}$ repeats with a period of $a-1$ (or one of its divisors) as $k$ varies, while every prime factor in Fibonacci-like recurrence relations repeats with a period of $a-1$, $a$, or $a+1$ (or divisors). Are there divisibility sequences where the factors have periods of e.g. $a-2$ or $a+2$? That is, anything other than $a$ or $a\pm 1$?
Let us consider a linear recurrence of order $d$, $$\tag1f(n)=a_1f(n-1)+\ldots + a_df(n-d) $$ where $a_1,\ldots, a_d\in \Bbb Z$ (with suitable start values $f(0),\ldots, f(d-1)$). As $d=1$ leads to constant $f$ or growth at least like $2^n$, we want to avoid that and assume henceforth that $d\ge 2$. Then the closed form is $$\tag2f(n):=c_1\lambda_1^n+\ldots+c_d\lambda_d^n,$$ where the $\lambda_\nu$ are the roots of the polynomial $X^d-a_1X^{d-1}-\ldots -a_d$ and the $c_\nu$ are determined by the inital values. (More precisely, if that polynomial has multiple roots, say $\lambda$ is a root of multiplicit $k$, then we have summands involving $\lambda^n, n\lambda^n,n^2\lambda^n,\ldots,n^{k-1}\lambda^n$ instead of $k$ instances of $\lambda^n$). It follows that $n\mapsto f(bn)$ has a similar closed form and thereby also follows a linear recurrence of order $d$.
In order to have a divisibilty sequence, we want $f(bn)$ to be a multiple of $f(b)$ for all $n$. In particular, $f(0)$ is a multiple of $f(b)$ for all $b$; if we want to exclude bounded $f$, this necessitates $$f(0)=0.$$ Then to have $f(b)\mid f(bn)$ for all $n$,
This is trivially the case if $d=2$.
What can we achieve with $d=2$?
First note that we can ignore the case of multiple roots when $d=2$: We have a double root $\lambda_1=\lambda_2=\frac{a_2}2$ iff $a_1^2=-4a_2$. If $c_1\lambda^n+c_2n\lambda^n$ is $0$ for $n=0$, we must have $c_1=0$, and in order to have unbounded $f$, we must have $a_1\ge 2$; we arrive either at a trivial solution of the form $f(n)=c_2n$ or at something that grows faster than the OP's $S_n$.
So we may assume different roots $\lambda_{1,2}=\frac{a_1\pm\sqrt{a_1^2+4a_2}}{2}$. If one of these is rational, then both are and in fact both are integers. This leads either to bounded $f$, or to $f$ growing at least like $2^n$. Hence $\lambda_{1,2}$ are irrational. From $f(0)=0$, we see $c_1=-c_2$ and conclude that the growth of $(2)$ is like $\max\{|\lambda_1|,|\lambda_2|\}^n$. But $\max\{|\lambda_1|,|\lambda_2|\}\ge \sqrt{|\lambda_1\lambda_2|}=\sqrt{|a_2|}$ and $\max\{|\lambda_1|,|\lambda_2|\}\ge \frac{|lambda_1+\lambda_2|}2=\frac{|a_1|}2$, so that in order to beat $\sqrt 2^n$, we need only consider $a_2=\pm1$ and $|a_1|\le 2$. The polynomials $X^2+2X+1$, $X^2+X+1$, $X^2+1$, $X^2-X+1$, $X^2-2X+1$, $X^2-1$ lead to roots of unity (sometimes double roots, which we excluded above), hence bounded $f$. The polynomials $X^2\pm 2X-1$ lead to growth like $(1+\sqrt 2)^n$, and the polynomials $X^2\pm X-1$ lead to Fibonacchi growth.
It follows that the OP's sequence $S_n$ is the slowest non-trivial order $2$ recurrence divisibility sequence. Upon closer inspection, $f(n)=0f(n-1)+2f(n-2)$ grows precisely as fast
What if $d>2$?
With larger $d$, we could clearly bring $\max\{|\lambda_\nu|\}$ closer to $1$. However, according to our necessary criterion above, we need to ensure that $f(2b)$ is a multiple of $f(b)$ for all $b$ (and more if $d>3$). In general, this seems not easily achievable. However, the solution obtained by stretching $\sqrt2^n$ apart, $$f(n)=0f(n-1)+0f(n-2)+2f(n-3) $$ of course grows only like $\sqrt[3]2^n$, and similarly for higher orders.