A maximum likelihood estimator has the nice property that it is invariant under transformations. This means that if $\theta_{MLE}$ is the MLE for $\theta$, then for a function $g(\theta)$, $g(\theta_{MLE})$ is the MLE for $g(\theta)$.
This can be directly applied to your problem. Hint: what is the MLE for $(\beta, \sigma)$?
The MLE is obtained by maximising the log-likelihood function, so the first thing you will want to do is have a look at this function. Using the density function for the Student T-distribution with known degrees-of-freedom $k \in \mathbb{N}$, you can write the log-likelihood as:
$$\ell_\mathbf{x,y} (\beta_0, \beta_1) = - \sum_{i=1}^n \ln \Big( k + (y_i - \beta_0 - \beta_1 x_i)^2 \Big).$$
The residuals $r_i = y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i$ under the MLE minimise $\sum_{i=1}^n \ln ( k + r_i^2 )$. As $k \rightarrow \infty$ we have $\ln ( k + r_i^2 ) = \ln(k) + \ln(1+r_i^2/k) \approx \ln(k) + r_i^2/k$ so that the residuals minimise $\sum_{i=1}^n r_i^2$ in the limit, which is the standard OLS solutions for normally distributed errors. Use of the T-distribution effectively dampens the effect of large residuals, through the above logarithmic transformation, and so the MLE is more tolerant of having some large residuals than in the normal case.
Finding the MLE: The MLE can be obtained via numerical maximisation of the log-likelihood using ordinary calculus techniques. The gradient of the log-likelihood is given by the partial derivatives:
$$\begin{equation} \begin{aligned}
\frac{\partial \ell_\mathbf{x,y}}{\partial \beta_0}(\beta_0, \beta_1)
&= \sum_{i=1}^n \frac{2 (y_i - \beta_0 - \beta_1 x_i)}{k+(y_i - \beta_0 - \beta_1 x_i)^2}, \\[6pt]
\frac{\partial \ell_\mathbf{x,y}}{\partial \beta_1}(\beta_0, \beta_1)
&= \sum_{i=1}^n \frac{2 x_i (y_i - \beta_0 - \beta_1 x_i)}{k+(y_i - \beta_0 - \beta_1 x_i)^2}. \\[6pt]
\end{aligned} \end{equation}$$
This leads to the score equations:
$$\begin{equation} \begin{aligned}
0 &= \sum_{i=1}^n \frac{2 (y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i)}{k+(y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i)^2}, \\[6pt]
0 &= \sum_{i=1}^n \frac{2 x_i (y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i)}{k+(y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i)^2}. \\[6pt]
\end{aligned} \end{equation}$$
These equations can be solved numerically using iterative techniques such as Newton-Raphson or gradient-descent methods, or other more complicated methods. The exact distribution of the MLE will be complicated, but for large samples the MLE should be normally distributed, according to standard large-sample theory.
Best Answer
In maximum likelihood estimation, we calculate
$$\hat \beta_{ML}: \sum \frac {\partial \ln f(\epsilon_i)}{\partial \beta} = \mathbf 0 \implies \sum \frac {f'(\epsilon_i)}{f(\epsilon_i)}\mathbf x_i = \mathbf 0$$
the last relation taking into account the linearity structure of the regression equation.
In comparison , the OLS estimator satisfies
$$\sum \epsilon_i\mathbf x_i = \mathbf 0$$
In order to obtain identical algebraic expressions for the slope coefficients we need to have a density for the error term such that
$$\frac {f'(\epsilon_i)}{f(\epsilon_i)} = \pm \;c\epsilon_i \implies f'(\epsilon_i)= \pm \;c\epsilon_if(\epsilon_i)$$
These are differential equations of the form $y' = \pm\; xy$ that have solutions
$$\int \frac 1 {y}dy = \pm \int x dx\implies \ln y = \pm\;\frac 12 x^2$$
$$ \implies y = f(\epsilon) = \exp\left \{\pm\;\frac 12 c\epsilon^2\right\}$$
Any function that has this kernel and integrates to unity over an appropriate domain, will make the MLE and OLS for the slope coefficients identical. Namely we are looking for
$$g(x)= A\exp\left \{\pm\;\frac 12 cx^2\right\} : \int_a^b g(x)dx =1$$
Is there such a $g$ that is not the normal density (or the half-normal or the derivative of the error function)?
Certainly. But one more thing one has to consider is the following: if one uses the plus sign in the exponent, and a symmetric support around zero for example, one will get a density that has a unique minimum in the middle, and two local maxima at the boundaries of the support.