Square-root matrix in Bimetric Gravity The Hasan-Rosen formulation of bimetric gravity can schematically be written as[1]:
\begin{equation}
\mathcal{L}_{bi} = \mathcal{L}_g + \mathcal{L}_f + \mathcal{L}_{int}
\end{equation}
where $\mathcal{L}_g$ is the Einstein-Hilbert action for metric $g_{\mu\nu}$ and similarly for metric $f_{\mu\nu}$. The interesting part really lies in the interaction term. This interaction term can schematically be written as:
$$ \mathcal{L}_{int} = \Sigma_{i=0}^{4}\  c_i \ e_i(\mathbb{X})  $$
Here: $e_i$ are the elementary symmetric polynomials. But the crux of this framework lies in the matrix $\mathbb{X}$, which is:
$$\mathbb{X} = \sqrt{g^{-1}f}$$
Question : Why the square-root? 
I understand that the simplest choice for interaction between two metrics is to contract them completely. Fine. But why is it put under a square-root?
[1]-https://arxiv.org/abs/1109.3515 (and there are many more references, but they all start with this construction)
 A: A partial answer to the question "why the square root" comes from realizing that the interaction terms are much simpler to write down in terms of vielbeins (aka tetrads/vierbeins/frame fields).  These vielbeins are roughly like a square root of the metric, and hence it is less surprising that the square root of $g^{-1} f$ appears when writing interaction terms between vielbeins.  
These statements can be made more precise.  The vielbeins are a collection of one-forms $e_\mu^a$ indexed by the internal Lorentz index $a$, and $\mu$ is just the spacetime index.  They are required to define an orthonormal frame,
$$g^{\mu\nu}e_\mu^a e_\nu^b = \eta^{ab},$$
with $\eta^{ab} = \text{diag}(-1,+1,+1,+1)$.  This equation also implies that 
$$e_\mu^a e_\nu^b \eta_{ab} = g_{\mu\nu},$$
so that $e_\mu^a$ acts like a square root of the metric, since a product of two of them gives back $g_{\mu\nu}$ (albeit with a necessary contraction with $\eta_{ab}$).  
To write the bimetric interactions, we also introduce an independent set of vielbeins $h_\mu^a$ for the second metric, satisfying
$$h_\mu^a h_\nu^b \eta_{ab} = f_{\mu\nu}$$
The bimetric interactions are then just given by
\begin{align}
L_1 &=e^a\wedge e^b\wedge e^c\wedge h^d \varepsilon_{abcd}\\
L_2 &=e^a\wedge e^b\wedge h^c\wedge h^d \varepsilon_{abcd}\\
L_3 &=e^a\wedge h^b\wedge h^c\wedge h^d \varepsilon_{abcd}
\end{align}
Here, $\varepsilon_{abcd}$ is simply the antisymmetric Levi-Civita symbol.  The spacetime indices are suppressed, and the wedge product refers to these spacetime indices.  We can convert these to Lagrangian densities by contracting with $\varepsilon^{\mu\nu\rho\sigma}$, the antisymmetric tensor density  (i.e. not divided by $\sqrt{-g}$).  Focusing just on $L_1$ for now, with a little work you can show that 
$$\frac1{4!}\varepsilon^{\mu\nu\rho\sigma}L^1_{\mu\nu\rho\sigma} = 3! \sqrt{-g} (e^{-1})_e^\alpha h_\alpha^e = 3! \sqrt{-g} \,\text{Tr}(e^{-1} h)$$ 
where we have introduced the inverse vielbein $(e^{-1})_e^\alpha = \eta_{ea}g^{\alpha\mu}e_\mu^a\equiv e^\alpha_e$, which acts like a square root of the inverse metric, $g^{\mu\nu} = e_a^\mu e_b^\nu \eta^{ab}$.  The other interactions involve symmetric polynomials of the matrix $e^{-1}h$.  So in particular, when written in terms of vielbeins, no square roots appear.  
The last thing we should ask is if $e^{-1} h$ really is the square root of $g^{-1} f$.  We can check this by squaring $e^{-1} h$, 
\begin{align}
e^{-1} h e^{-1} h &=e^\mu_a h^a_\alpha e^\alpha_bh^b_\nu \\
&= g^{\mu\beta}e_\beta^c\eta_{ca}h^a_\alpha e^\alpha_b h^b_\nu
\end{align}
At this point, we would like to commute the $e^c_\beta$ past the $h^a_\alpha$, in order to contract it with another $e$ vielbein.  It turns out that we can do this if we assume what is known as the symmetric vielbein condition,
$$\eta_{ca} e^a_\beta h^c_\alpha = \eta_{ca}e^a_\alpha h^c_\beta.$$
Using this, we can complete the calculation,
\begin{align}
e^{-1} h e^{-1} h &= g^{\mu\beta} h^a_\beta \eta_{ac} e^c_\alpha e^\alpha_b h^b_\nu \\
&= g^{\mu\beta}h^a_\beta \eta_{ab} h^b_\nu\\
&= g^{\mu\beta} f_{\beta\nu} \;\checkmark
\end{align}
Finally, you could ask whether imposing the symmetric vielbein condition is justified.  This issue is discussed in this paper https://arxiv.org/abs/1208.4493, where the conclusion seems to be that the condition is justified precisely when $g^{-1} f$ has real square roots.  
