I have the following optimization problem.
\begin{align} \text{Maximize } & \sum_{k=1}^P x_{i,k} \mu_{k} \\ \text{Subject to } & \sum_{k=1}^P\sum_{l=1}^P x_{i,k}x_{i,l} \sigma_{k, l} \geq \epsilon_i \\ & \sum_{k=1}^P\sum_{l=1}^P x_{i,k}x_{j,l}^*\sigma_{k, l} \leq \gamma_{i, j} \; \forall \; 1 \leq j \leq i - 1 \\ & \sum_{k=1}^P x_{i, k} = K \end{align}
With decision variables $x_{i, k} \in \{0, 1\} \; \forall k \in \{1, 2, \dots, P\}$ and all other variables constant.
My problem is that the first constraint is not linear. In the paper which this program is adapted from (page 12, eq 2.14), they state "We can introduce auxiliary variables to make this formulation linear as is standard in many integer programming formulations.", but I don't know enough about integer linear programming to know how to do that.
How do I make this constraint linear?