Tukukkk

If $a_n$ satisfies the linear recurrence relation $a_n = sum_{i=1}^k c_i a_{n-i}$ for some constants $c_i$, then is there an easy way to find a linear recurrence relation for $b_n = a_n^2$ ?

For example, if $a_n = a_{n-1} + a_{n-3}$, then $b_n=a_n^2$ seems to satisfy $b_n=b_{n-1}+b_{n-2}+3b_{n-3}+b_{n-4}-b_{n-5}-b_{n-6}$.

Yes. Take the companion matrix $M$ of the characteristic polynomial of your original recurrence. Then the squared recurrence satisfies a recurrence with characteristic polynomial the characteristic polynomial of the symmetric square $S^2(M)$ of $M$. If the original recurrence has order $n$ this new recurrence has order ${n+1 choose 2}$, and its coefficients are polynomials (depending only on $n$) in the coefficients of the old recurrence.

To prove this it suffices to consider the case where the characteristic polynomial has distinct roots, by a density argument. Then $a_n$ is a linear combination of the sequences $lambda_i^n$ where $lambda_i$ are the roots of the characteristic polynomial. So $a_n^2$ is a linear combination of the sequences $(lambda_i lambda_j)^n$ (and we may have $i = j$), and $S^2(M)$ has the $lambda_i lambda_j$ as its eigenvalues. In general these are all also distinct, so the characteristic polynomial of $S^2(M)$ is minimal with this property and it's not possible to reduce the order further than this in general.

The same argument shows that for $k^{th}$ powers we can use the symmetric powers $S^k(M)$, and that for general products $a_n b_n$ of sequences satisfying linear recurrences we can use tensor / Kronecker products of the companion matrices of their characteristic polynomials.

T. Brown and P.J. Shiue's paper here might be of interest as a first reference. In the introduction they mention that if $a_n$ is a second-order sequence then the sequence
of squares $a_{n}^2$ is a third-order sequence. They go on to show necessary conditions for the squares sequence to be a second-order sequence when $a_n$ is a homogeneous sequence.

In the paper of Cooper and Kennedy here, (section 5) they give an order six linear recurrence relation for the square of a third order linear recurrence relation (as appears in your example):

$$x^2_n = (a^2 + b)x^2_{n−1} + (a^2b + b^2 + ac)x^2_{n−2} + (a^3c + 4abc − b^3 +2c^2)x^2_{n−3}+(−ab^2c + a^2c^2 − bc^2)x^2_{n−4} + (b^2c^2 − ac^3)x^2_{n−5} − c^4x^2_{n−6}$$
where $x_n = ax_{n−1} + bx_{n−2} + cx_{n−3}$.
For similar questions where squares have been replaced with higher powers, the paper here by Stinchcombe might be interesting.

Edit: Qiaochu Yuan has provided the correct order for squares. Higher powers are addressed using the same argument in Theorem 3 of Stinchcombe's paper above. For convenience it is stated here:

Question: What is the order of the recurrence relation $y_{n} =x^l_{n}$, for $x_{n}$ a recurrence relation of order $k$?
The degree of the corresponding characteristic polynomial is counted by the
number of elements in $B_{l}$, where
$$ B_{l} = {(i_{1},...,i_{k}) | text{ each }i_j text{ is a nonnegative integer and } i_1+...+i_{k}=l }.$$
Given a value of $k$, define $S(k, l) =|B_{l}|$, then

Theorem 3: $S(k, l)$ obeys the relations: $S(k, l) = k$ for all $k$, $S(1, l) = l$ for all $l$, and $S(k,l) =S(k-1,l) + S(k, l-1)$ for every $k$ and $l$.