Eigen equivalent to arithmetic with Numpy Brodcasting (3D)

In numpy i have the following

X1 = np.arange(1, 6).reshape(-1, 1)
X1 = np.hstack([X1, X1])

# Output: X1
# array([ [1, 1],
#         [2, 2],
#         [3, 3],
#         [4, 4],
#         [5, 5] ])

X2 = X1[:, np.newaxis, :]
X3 = X1[np.newaxis, :, :]

# Output: Shapes
# X2.shape = (5, 1, 2)
# X3.shape = (1, 5, 2)

X4 = (X2 - X3) ** 2

# Output: X4 
# X4.shape = (5, 5, 2)
# X4[:, :, 0] = 
# array([ [ 0,  1,  4,  9, 16],
#         [ 1,  0,  1,  4,  9],
#         [ 4,  1,  0,  1,  4],
#         [ 9,  4,  1,  0,  1],
#         [16,  9,  4,  1,  0] ], dtype=int32)
# X4[:, :, 1] = 
# array([ [ 0,  1,  4,  9, 16],
#         [ 1,  0,  1,  4,  9],
#         [ 4,  1,  0,  1,  4],
#         [ 9,  4,  1,  0,  1],
#         [16,  9,  4,  1,  0] ], dtype=int32)

I would need this because I will later have to divide X4 by another (Row) vector with the same number of columns as X1. I understand Eigen has its broadcasting operations (https://eigen.tuxfamily.org/dox/group__TutorialReductionsVisitorsBroadcasting.html) however my numpy implementation works with 3D matrices.

My solution was to use for loops, however is there a more efficient way to do this?

typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Matrix;
void test(const Matrix& X1, const Matrix& X2, std::vector<Matrix>& D) {

    for (int i = 0; i < X1.cols(); i++) {
        Matrix tmp(X1.rows(), X2.rows());
        for (int j = 0; j < X1.rows(); j++) {
            for (int k = 0; k < X2.rows(); k++) {
                tmp(j, k) = X1.col(i)(j) - X2.col(i)(k);
            }
        }
        D.push_back(pow(tmp.array(), 2));
    }
}