Let $X^a$ and $X^b$ be geometric Brownian motions, i.e. for any $i\in \{a,b\}$, $$ dX_t^i/X_t^i=\mu_idt+\sigma_idB^i_t $$ with $X_0^i=x_i>0$, where $B^a$ and $B^b$ are independent Brownian motions. It is clear that $X_t^i=x_ie^{(\mu_i-\sigma_i^2/2)t+\sigma_iB^i_t}$ a.s. and $\mathbb{E}[X_t^i]=x_ie^{\mu_i t}$.
Question: Let $p\in(0,1)$, calculate $\mathbb{E}[(X_t^a+X_t^b)^p]$.
My attempt: \begin{align*} \mathbb{E}[(X_t^a+X_t^b)^p]&=\mathbb{E}\left[(x_ae^{(\mu_a-\sigma_a^2/2)t+\sigma_aB^a_t}+x_be^{(\mu_b-\sigma_b^2/2)t+\sigma_bB^b_t})^p\right]\\ &=\mathbb{E}\left[(x_ae^{(\mu_a-\sigma_a^2/2)t+\sigma_a\sqrt{t}Z^a}+x_be^{(\mu_b-\sigma_b^2/2)t+\sigma_b\sqrt{t}Z^b})^p\right]\\ &=\mathbb{E}\left[\int_{\mathbb{R}^2}(x_ae^{(\mu_a-\sigma_a^2/2)t+\sigma_a\sqrt{t}u}+x_be^{(\mu_b-\sigma_b^2/2)t+\sigma_b\sqrt{t}v})^pf(du)f(dv)\right] \end{align*} where $Z^a$ and $Z^b$ are standard normal random variables and the function $f$ is the density of standard normal distribution. I tried to calculate this integral by tools such as Mathematica but none worked. Alternatively, one can try to solve the PDE from the martingale approach i.e. setting $E[(X_t^a+X_t^b)^p|\mathcal{F}_s]=:M_s$ and the PDE is yielded using Ito formula. However, PDE is hard to solve.