To find the weaks forms we use the identity
$$
\nabla\cdot(\rho\,\nabla\phi) = \rho\,\nabla^2\phi + \nabla\rho \cdot \nabla\phi
$$
and its vector counterpart
$$
\nabla\cdot(\varphi\,\mathbf{v}) = \varphi\,(\nabla\cdot\mathbf{v}) + (\nabla\varphi) \cdot \mathbf{v} \,.
$$
Then the first equation can be written as
$$
\partial_{t}\rho+2\,\nabla\cdot(\rho\,\nabla\phi)+2\sigma\rho\left(\rho-1\right)=0 \,.
$$
If we consider only the spatial terms and use a weighting function $w$, the weak form is
$$
\int_\Omega (\partial_{t}\rho)\,w\,d\Omega + 2\int_\Omega \left[\nabla\cdot(\rho\,\nabla\phi)+\sigma\rho\left(\rho-1\right)\right]w\,d\Omega=0 \,.
$$
Using the vector identity above, we get
$$
\int_\Omega (\partial_{t}\rho)\,w\,d\Omega + 2\int_\Omega \left[\nabla\cdot(w\,\rho\,\nabla\phi)-\nabla w\cdot(\rho\,\nabla\phi) + w\sigma\rho\left(\rho-1\right)\right]\,d\Omega=0 \,.
$$
Apply the divergence theorem to get
$$
\int_\Omega (\partial_{t}\rho)\,w\,d\Omega + 2\int_\Gamma (w\,\rho\,\nabla\phi)\cdot\mathbf{n}\,d\Gamma - 2\int_\Omega \left[\nabla w\cdot(\rho\,\nabla\phi) - w\sigma\rho\left(\rho-1\right)\right]\,d\Omega=0 \,.
$$
Using the same identity as before, we also have
$$
\nabla\phi\cdot\nabla\phi = \nabla\cdot(\phi\,\nabla\phi) - \phi\,\nabla\cdot(\nabla\phi)
$$
and
$$
\nabla\cdot\left(\frac{1}{\rho}\,\nabla\rho\right) = \frac{\nabla^2\rho}{\rho} -\frac{\nabla\rho \cdot \nabla\rho}{\rho^2}\,.
$$
Assuming that the extra factor of $1/2$ in the second equation is a typo, we can write that equation as
$$
\partial_{t}\phi+\nabla\cdot(\phi\,\nabla\phi) - \phi\,\nabla\cdot(\nabla\phi)-\frac{1}{2}\nabla\cdot\left(\frac{1}{\rho}\,\nabla\rho\right)+\rho-1=0
$$
It is here that the appropriate choice of terms to consider at the boundaries becomes important. The choice is typically determined by the ease of application of the boundary conditions and how the two sets of equations are coupled in the numerical method.