How do I convert this to a linear integer program?

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?

Let me introduce variable $y_{i,k,l}= x_{i,k}x_{i,l}$,

We want $$y_{i,k,l} = \begin{cases} 1, & x_{i,k}=1 \text{ and } x_{i,l}=1\\ 0, & \text{Otherwise} \end{cases}$$

Let's replace the first constraint by:

$$y_{i,k,l} \geq x_{i,k} + x_{i,l}-1$$

$$y_{i,k,l} \leq x_{i,k}$$

$$y_{i,k,l} \leq x_{i,l}$$

$$0 \leq y_{i,k,l} \leq 1$$

$$\sum_{k=1}^P\sum_{l=1}^P y_{i,k,l} \sigma_{k, l} \geq \epsilon_i$$