When does the sum of two $t$-distributed random variables follow a $t$ distribution? In the scope of a project, I need to find the sum of two independent $t$-distributions.
I know that in the general case, the sum of two $t$-distributed random variables is not $t$-distributed. However, I would like to know what might be the special cases where the sum of two $t$-distributions is a $t$ distribution.
It would also be nice if someone can give me an example of two such $t$-distributions.

Characteristic functions are handy for studying sums of independent variables, because the c.f. of the sum is the product of the c.f.s of the addends and c.f.s determine the distribution.
The c.f. of a Student t variable with parameter ("degrees of freedom") $\nu \ge 1$ is, up to a constant of proportionality $C(\nu)$ we don't need to consider, given by
$$\phi_\nu(t) = C(\nu) K_{\nu/2}(|t|\sqrt{\nu})\,|t|^{\nu/2}\tag{*}$$
where $K$ is the Bessel Function of the Second Kind.  I promise we won't need to learn much about this function in the following analysis.
Interpreting the question a little more broadly to accommodate a possible change of scale in the result, it asks for all solutions $\nu,\kappa,\mu,\sigma$ to the equation
$$\phi_\nu\,\phi_\kappa= \phi_\mu\circ \lambda_\sigma$$
where all the subscripts are Student t parameters (necessarily positive) and the new scale $\sigma$ is positive; $\lambda_\sigma (t) = \sigma t$ is the change of scale.  According to $(*),$ this equation is equivalent to
$$K_{\nu/2}(|t|\sqrt{\nu})\, K_{\kappa/2}(|t|\sqrt{\kappa}) |t|^{(\nu + \kappa)/2}\ \propto\ K_{\mu/2}(|t|\sigma \sqrt{\mu}) |t|^{\mu/2}.\tag{**}$$
As $|t|$ grows large, $K_\alpha(|t|)$ is asymptotically proportional to
$$K_\alpha(|t|) \ \sim\ C(\alpha)e^{-|t|}|t|^{-1/2}\left(1 +  \frac{4\alpha^2-1}{8|t|} + O\left(|t|^{-2}\right)\right).$$
Comparing the leading terms of $(**)$ shows us that
$$e^{-(\sqrt\nu + \sqrt\kappa)|t|}|t|^{(\nu + \kappa)/2-1}\ \propto\ e^{-\sigma\sqrt\mu|t|}|t|^{(\mu-1)/2},$$
equivalent to $\nu+\kappa-1=\mu$ and $\sqrt{\nu}+\sqrt{\kappa} = \sigma\sqrt{\mu}.$  Consequently, one or both of the original random variables must have a Student $t(1)$ distribution.  Moreover, taking $\kappa=1$ in that case, whence $\mu=\nu,$ we find $\sigma\sqrt{\nu} = \sqrt{\nu}+1,$ entailing
$$\sigma = 1 + \frac{1}{\sqrt{\nu}}.$$
Comparing the $|t|^{-1}$ terms in the asymptotic expansion shows
$$\frac{4\nu^2-1}{8|t|} =\frac{4\nu^2-1}{8\sigma|t|} + O\left(|t|^{-2}\right)$$
which is possible only when $1 = \sigma = 1 + 1/\sqrt{\nu}$ or $4\nu^2-1=0:$ the only solution is $\nu=1.$

The sum of two independent Student t variables has a Student t distribution (up to scale) only when both variables have one degree of freedom; and in that case, the resulting distribution has one degree of freedom and a scale factor of $2.$

As a check, here is a QQ plot of $|X+Y|$ for a million independent realizations of a Student $t(2)$ variable $Y$ and Student $t(1)$ variable $X$ (shown on a log-log scale due to the heavy tail).

(To avoid plotting all million values, only the 0.5%, 1%, ..., 99.5% quantiles of the data are shown.)
The agreement with a $\sigma$ - scaled Student $t(2)$ distribution is good in the middle, but grows increasingly worse out in the tail: adding the heavy-tailed $X$ to $Y$ has given the sum a heavier tail, even after scaling.  It is this deviation at the upper right that the analysis of the $O(|t|^{-1})$ term in the c.f. revealed.
