Optimize the product of three binary variables: convex relaxation and integer solutions

Question

I am working on a problem related to graphs and I formulated the problem as follows:

$$\max_{y,e} \sum_{i=1}^n (c_i^1 y_i^1 + c_i^2 y_i^2) + \sum_{(i,j)\in E}(d^1y_i^1y_j^1e_{ij} +d^2y_i^2y_j^2e_{ij})-\sum_{(i,j)\in E}a_{ij}e_{ij}$$ s.t. $$y_i^1 + y_i^2 = 1,\ i=1,2,\cdots,n$$ $$\sum_{(i,j)\in E} e_{ij} \geq D$$

The variables are $(\cdots, y_i^1,y_i^2,\cdots)$ for $i=1,2,\cdots,n$ and $e_{ij}$ for $(i,j)$ in an edge set $E$. All variables are binary.

$(c_i^1,c_i^2)$ are (unrestricted) constants and $d^1,d^2,a_{ij},D$ are non-negative constants.

The above objective function is nonlinear and the variables are binary. So I used the standard method to linearize the objective:

Let $z_{ij}^1 = y_i^1y_j^1e_{ij}$ and $z_{ij}^2 = y_i^2y_j^2e_{ij}$ and add constraints $z_{ij}^1 \leq y_i^1,z_{ij}^1 \leq y_j^1,z_{ij}^1 \leq e_{ij},z_{ij}^1 \geq y_i^1 + y_j^1 + e_{ij} -2$ and similarly for each $z_{ij}^2$.

Second, I relax the integer constraints and let every variable be within $[0,1]$.

The resulting problem is a linear program, which can be solved optimally. However, the solutions are fractional.

Below are some issues.

1. Is there any condition (for the form of the objective function) with which the linear program is guaranteed to produce integer solutions? I met with similar problems where the nonlinear terms in the objective only involve products of two binary variables and using the same linearization technique, the solutions are integral. Also, I removed the terms $-\sum_{(i,j)\in E} a_{ij}e_{ij}$ in the original objective functions and then it seems that the solutions are all integral.

2. Is there any other technique to linearize the product of three binary variables? What are the standard method(s) to deal with such optimization problems?

3. Since the solutions are fractional, what are the standard method(s) to round the results to get integer solutions? And are there any approximation guarantees for such rounding methods? I tried a few methods, such as setting a threshold (fractional values above the threshold are set to 1) and ranking the fractional solutions (pick the top-K to be 1).