Let $Y_1\sim SN(\mu_1,\sigma_1^2,\lambda)$ and $Y_2\sim N(\mu_2,\sigma_2^2)$ independents. Show that $Y_1+Y_2$ have a skew-normal distribution and find the parameters of this distribution.
Since the random variables are independent I tried to use convolution. Let $Z=Y_1+Y_2$
$$f_Z(z)=\int_{-\infty}^{\infty}2\phi(y_1|\mu_1,\sigma_1)\Phi\Big(\lambda(\frac{y_1-\mu_1}{\sigma_1})\Big)\phi(z-y_1|\mu_2,\sigma_2^2)\,\text{d}y_1$$
Here $\phi()$ and $\Phi()$ are the standard normal pdf and cdf, respectively.
$$f_Z(z)=\int_{-\infty}^{\infty}2\frac{1}{\sqrt{2\pi\sigma_1}}\frac{1}{\sqrt{2\pi\sigma_2}}exp\Big(-\frac{1}{2\sigma_1^2}(y_1-\mu)^2-\frac{1}{2\sigma_2^2}((z-y_1)^2-\mu)^2\Big)\Phi\Big(\lambda(\frac{y_1-\mu_1}{\sigma_1})\Big)\,\text{d}y_1$$
For simplified notations, let $k=2\frac{1}{\sqrt{2\pi\sigma_1}}\frac{1}{\sqrt{2\pi\sigma_2}}$
\begin{align*}f_Z(z)&=k\int_{-\infty}^{\infty}\exp\Big(\frac{-1}{2\sigma_1^2\sigma_2^2}\Big(\sigma_1^2(y_1-\mu_1)^2+\sigma_2^2((z-y_1)-\mu_2)^2\Big)\Big)\Phi\Big(\lambda(\frac{y_1-\mu_1}{\sigma_1})\Big)\,\text{d}y_1\\ &=k\int_{-\infty}^{\infty}\exp\Big(\frac{-1}{2\sigma_1^2\sigma_2^2}\Big(\sigma_2^2(y_1^2-2y_1\mu_1+\mu_1)+\sigma_1^2((z-y_1)^2-2(z-y_1)\mu_2+\mu_2^2)\Big)\Big)\\&\quad\times\Phi\Big(\lambda(\frac{y_1-\mu_1}{\sigma_1})\Big)\,\text{d}y_1=k\int_{-\infty}^{\infty} \exp\\&\Big(\frac{-1}{2\sigma_1^2\sigma_2^2}\Big(\sigma_2^2(y_1^2-2y_1\mu_1+\mu_1)+\sigma_1^2(z^2-2zy_1+y_1^2-2z\mu_2+2y_1\mu_2+\mu_2^2)\Big)\Big)\\&\quad\times\Phi\Big(\lambda(\frac{y_1-\mu_1}{\sigma_1})\Big)\,\text{d}y_1 \end{align*}
But I'm stuck at this point.
EDIT: Following the suggestions in the comments, taking $\mu_1=\mu_2=0$ and $\sigma_1^2=\sigma_2^2=1$ \begin{align*} &\int_{-\infty}^\infty 2\frac{1}{\sqrt{2\pi}}\frac{1}{\sqrt{2\pi}}\exp\Big(-\frac{1}{2}[y_1^2+z^2-2zy_1+y_1^2]\Big)\Phi(\lambda y_1)dy_1 \\&\int_{-\infty}^\infty 2\frac{1}{\sqrt{2\pi}}\frac{1}{\sqrt{2\pi}}\exp\Big(-\frac{1}{2}y_1^2\Big)\Phi(\lambda y_1) \exp\Big(-\frac{1}{2}(z-y_1)^2\Big)dy_1\end{align*}
is skew-normal.