Given an $ntimes ntimes n$ cube, what is the largest number of $1times 1times 1$ blocks that a plane can cut through?

Let's say the origin of the coordinate system is the center of $n times n times n$ cube.
Also, each side of the cube is parallel to each axis of the coordinate system.
Let the plane $P$ cuts the cube. equation of plane $P$ is
$$ P : ax+by+cz+d=0 ; ; (a geq 0 , bgeq 0 , c geq 0 , a^2+b^2+c^2=1)$$
(Because it's the same thing when you rotate it, $(a geq 0 , bgeq 0 , c geq 0)$)

Case 1: $n$ is odd
the area of each block $B_{k m l}$ is given as follows.
(For three integers $k,m,l$ that satisfy $|k|,|m|,|l| leq left [ frac{n}{2} right ] $ ) $$B_{k m l} := left { (x,y,z) : |x-k| leq 1/2 , |y-m| leq 1/2 , |z-l| leq 1/2 right }$$
In thin case, Set $E_{k m l}$ , consists of all vertices in block $B_{k m l}$ is given as follows
$$E_{k m l} := left { (x,y,z) : x= kpm1/2 , y=mpm1/2 , z = l pm 1/2 right } $$

If block $B_{k m l}$ is cut by plane $P$, some $ mathbf{a},mathbf{b} in E_{k m l}$,$f(mathbf{a})f(mathbf{b})<0$ where $f(x,y,z) := ax +by +cz+d$
For $mathbf{a} in E_{k m l}$ , define the minimum value of $f(mathbf{a})$ as $min(k,m,l)$ , and define the maximum value of $f(mathbf{a})$ as $max(k,m,l)$.
Then, $$ min(k,m,l) = f(k,m,l) -frac{1}{2}(a+b+c) ; , ; max(k,m,l) = f(k,m,l) + frac{1}{2}(a+b+c) $$
Therefore, $ |f(k,m,l)| < frac{1}{2}(a+b+c)$ is a necessary and sufficient condition for block $B_{k m l}$ to be cut by plane $P$.
This condition means, point $(k,m,l)$ must be between planes $P^{+} : ax+by+cz+d = frac{1}{2}(a+b+c)$ and $P^{-} : ax+by+cz+d = -frac{1}{2}(a+b+c)$.
Note that The distance between the two planes is $sqrt{3}$.

I guess it should be $d=0$ , And I think there will be a value of appropriate $(a,b,c)$ (regardless of the value of $n$).
I'm sorry but I don't know the specific way to prove this.

The above explanation are too complicated for me. I made an error in counting so here is how you do it and from this diagram some generalization may be made. The picture is a top view of the 3X3X3 cube. The diagonal lines are the cut intersections with the borders of the layers of 3X3 cubes. The numbers represents the cubes that are cut in each layer - 1 for lower layer, 2 for middle layer, and 3 for top layer.

Bottom (1) and top (3) layers have 6 cubes cut and the middle (2) layer have 7 - total 19. I do not see a way to make 20.

The hexagonal cross-section midway between diagonally ooposite vertices has side-length $n/sqrt2$ and area $(3sqrt3/4)n^2$. It cuts a cube when the cube's centre is within $sqrt3/2$ of the plane. The available volume is then $ (9/4)n^2$, so the leading order number of cut cubes might be $(9/4)n^2$.
Let a normal to the plane be $(a,b,c)$. By symmetry, we can assume that $a,b,c$ are all positive.
Whenever the cross-section is a hexagon, the normal vector $(a,b,c)$ satisfies the triangle inequalities $$alt b+c \blt a+c \c lt a+b$$
So we can write $a=u+v, b=u+w, c=v+w$ for positive $u,v,w$.
The cross-sectional area for a given normal is greatest when the plane goes through the centre of the $n×n×n$ cube.
The available volume for the cut cubes' centres is
$$left(frac{(a+b+c)^3-2(a+b+c)(a^2+b^2+c^2)}{4abc}right)n^2 \ =left(frac{2(uv+uw+vw)(u+v+w)}{(u+v)(u+w)(v+w)}right)n^2 \ =left(frac94-frac{u(v-w)^2+v(u-w)^2+w(u-v)^2}{4(u+v)(u+w)(v+w)}right)n^2$$ So the available volume is at most $9n^2/4$ when the cross-section is a hexagon.

This answer solves the $3times 3times 3$ case and makes a conjecture about further cases.

To give an answer, first imagine how we can create the given $ntimes ntimes n$ cube in the first place: take all of $mathbb R^3$. Draw $(n+1)$ equally spaced parallel planes. Discard everything lying "outside" of the outermost two planes and imagine the space within to be cut by each of the remaining $(n-1)$ planes. Repeat this process for a set of planes perpendicular to the original set, and then for a set of planes perpendicular to both sets so far each spaced equally to the first.

Note that "perpendicular" is irrelevant here as is the equality of the spacing between the three sets, since the problem only references linear structure - so as long as we picked the orientations of the planes to be independent and keep the spacing within each set constant, the problem is unchanged.

The trick is to choose whichever plane we wish to use to slice the cube first and then to perform the above procedure and see what happens to the plane. In particular, after the first two sets of slices, the plane will have been reduced and cut into an $ntimes n$ grid of parallelograms - and, again, only linear structure being relevant, we may as well reduce to the following question:

Suppose we have an $ntimes n$ grid of squares. Draw a set of $(n+1)$ parallel and equally spaced lines. Discard all squares fully outside of the bounds of these lines and imagine the squares to be cut at each of the remaining $(n-1)$ lines. How many regions could there remain?

This question seems more approachable - as it happens on a 2D grid rather than a 3D space. A lot of subtleties happen when one tries to solve the above question, though - you should certainly not have any of the additional lines pass through any corners of squares, as perturbing any cut having this property would yield more pieces. Moreover, you can express the number of pieces cut as "the number of squares not entirely discarded plus the number of line segments cut off from the intermediate lines by the squares."

You certainly can't do better than cutting $n^2 + (n-1)(2n-1)=3n^2-3n+1$ regions by following out the logic above, but achieving this would require that no square is entirely discarded, but each intermediate cut cuts through the maximum of $(2n-1)$ interior squares - which is clearly impossible for large $n$.

I might guess that the optimal configuration for $ngeq 3$ is to take the longest diagonal on the $ntimes n$ squares and draw the further $(n+1)$ lines to all hit every point on that diagonal and so that this each square within two squares of the diagonal has at least part of itself between the outer bounding lines and so that every intermediate line hits every square on the diagonal without being the exact diagonal - meaning each intermediate line crosses $(2n-1)$ squares and that $n+2(n-1)+2(n-2)$ squares are not entirely discarded and $(n-1)(2n-1)$ are cut by the intermediate lines - for a total of $2n^2+2n - 5$ regions left - i.e. that a plane across a $ntimes n times n$ cube can hit at least $2n^2+2n-5$ of the $1times 1times 1$ cubes. This might be optimal, but it's unclear if widening the distance between the outer lines to include more squares at least partially might offset that some lines would then create fewer new regions - and the reasoning to figure that out seems really sensitive, since no matter what you do, it seems like, no matter what you do, you'll remain on the order of $2n^2$ with only lower order terms up for grabs.

Notice that the lower bound and upper bound are both equal to $19$ when $n=3$ - so this is the answer for a $3times 3times 3$ cube and a conjecture for larger cubes. For concreteness, if we assume that this is the cube $[-3,3]times [0,3]times [0,3]$, a plane achieving this maximum is defined by $z = x+y-frac{3}2$, noting every relevant square of the $x$-$y$ plane lies at least partially in the region $0leq x+y - frac{3}2leq 3$ - so a cube in every $z$ column is included - and the lines $x+y-frac{3}2=1$ and $x+y-frac{3}2=2$ each hit five squares - contributing an extra cube for each of these incidences, for a total of $10$ cubes (or, specifically: two corner columns have $1$ cube hit each, the four middle-of-edge columns have $2$ cubes hit each, and the three diagonal columns each get $3$ cubes hit each - for a total of $19$ cubes hit by the plane).

---

Edit: Some computational results: if we consider only the planes of the form $x+y+kcdot z = (k+2)n/2$ - which are planes passing through the center pivoting around a certain axis (chosen so that in the diagram of squares, the added lines are diagonal - though there's no formal reason to believe this is optimal) - we can actually just use a computer to check what the optimal $k$ are. The suggested optimal setup above is not optimal for all $n$ (and neither is the suggestion of choosing $k=1$).

For $n=3$, a maximum of $19$ cubes hit by such planes is achieved for $2/3 < k < 2$. For $n=4$, a maximum of $35$ cubes can be hit for $1/2 < k < 1$. For $n=5$ a maximum of $57$ cubes can be hit for $5/4 < k < 4/3$. For $n=6$ a maximum of $81$ are hit for $2/3 < k < 1$. For $n=7$ a maximum of $113$ cubes can be hit for $8/7 < k < 5/4$. For $n=8$ we get a maximum of $145$ for $3/4 < k < 1$. For $n=9$, we get a maximum of $187$ cubes hit for $10/9 < k < 9/8$. There seem to be some patterns, but the plots of the number of cubes hits vs. the slope are very uneven, jumping up and down seemingly randomly and clearly depending on parity. This problem might not be as clear cut as I had thought - no idea how to solve it in general.

Given a cube $n times n times n$ or $[0,, n]^3$ we want to find the plane $ax+by+cz=d$ which crosses the highest number of unitary cubes inside $[0,, n]^3$, and find that number.

We individuate a single unit cube as $[x_k,, x_k+1] times [y_j,, y_j+1] times [z_l,, z_l+1]$, with $j,k,l in [0, , n-1]$.

The cubes crossed by the plane will be those for which $$ eqalign{ & ax_{,k} + by_{,j} + cz_{,l} < d < aleft( {x_{,k} + 1} right) + bleft( {y_{,j} + 1} right) + cleft( {z_{,l} + 1} right)quad Rightarrow cr & Rightarrow quad d - left( {a + b + c} right) < ax_{,k} + by_{,j} + cz_{,l} < dquad Rightarrow cr & Rightarrow quad {d over {a + b + c}} - 1 < {{ax_{,k} + by_{,j} + cz_{,l} } over {a + b + c}} < {d over {a + b + c}} cr} $$

Consider $x_k$ as the realization a uniform discrete random variable $x$ on the support $[0,, n-1]$, with probability $1/n$, mean $(n-1)/2$ and variance $(n^2-1)/12$ .
Same for $y, , z$.

Their weighted sum $$ {{ax_{,k} + by_{,j} + cz_{,l} } over {a + b + c}} $$ will have mean, mode and median at $(n-1)/2$ and variance $$ sigma ^2 = {{a^2 + b^2 + c^2 } over {left( {a + b + c} right)^2 }}left( {{{n^2 - 1} over {12}}} right) $$

Clearly the less is the variance the larger is the portion of the pmf satisfying the inequality given above, since the inequality's gauge is constant at $1$.
And the variance is clearly minimum for equal weights.

So we arrive to consider the inequality $$ bbox[lightyellow] { left{ matrix{ x_{,k} ,y_{,j} ,z_{,l} ,n,s in mathbb Z hfill cr d in mathbb R hfill cr 0 le x_{,k} ,y_{,j} ,z_{,l} le n - 1 hfill cr d - 3 < x_{,k} + y_{,j} + z_{,l} = s < d hfill cr} right. tag{1}}$$

Now, the number of points on the diagonal plane of a $m$-D cube $$N_{,b} (s,r,m) = text{No}text{. of solutions to};left{ begin{gathered} 0 leqslant text{integer }x_{,j} leqslant r hfill \ x_{,1} + x_{,2} + cdots + x_{,m} = s hfill \ end{gathered} right.$$ is given by $$ bbox[lightyellow] { N_b (s,r,m)quad left| {;0 leqslant text{integers }s,m,r} right.quad = sumlimits_{left( {0, leqslant } right),,k,,left( { leqslant ,frac{s}{r+1}, leqslant ,m} right)} {left( { - 1} right)^k binom{m}{k} binom { s + m - 1 - kleft( {r + 1} right) } { s - kleft( {r + 1} right)} } tag{2.a}}$$ as explained in this post.

Moreover the number of points on or below the diagonal plane are $$ bbox[lightyellow] { eqalign{ & M_b (s,r,m)quad left| {;0 le {rm integers }s,m,r} right.quad = cr & = sumlimits_{0, le ,,j,, le ,s} {N_b (s,r,m)} = cr & = sumlimits_{left( {0, le } right),,k,,left( { le ,{s over {r + 1}}, le ,m} right)} {left( { - 1} right)^k left( matrix{ m hfill cr k hfill cr} right)left( matrix{ s + m - kleft( {r + 1} right) cr s - kleft( {r + 1} right) cr} right)} cr} tag{2.b}}$$

At this point we need the help of a graphical visualization to grasp the behaviour of the inequality 1) wrt $N_b$

The sketch represents the histograms of $N_{,b} (s,n-1,3)$ for $n=3$ and $n=4$.
$N_{,b} (s,n-1,3)/n^3$ is the pmf of the sum $s$ of the three uniform discrete random variables.
The sketch shows that the maximum portion of the histogram is intercepted when the the gauge of width $3$ of the inequality is almost centered around the mean.
That is actually so when n is odd, while for even $n$ we shall shift the gauge slightly to the left (or to the right).
Alas, the formula for $N_b$ is only valid for integral parameters (rewriting the binomial through gamma produces a discontinuous function).

We can circumvent the above and uniform the inequality by introducing a fixed $1/2$ shift from the mean and then rewriting the inequality as $$ eqalign{ & d - 3 < x_{,k} + y_{,j} + z_{,l} = s < dquad Rightarrow cr & Rightarrow quad 3{{n - 1} over 2} - 3/2 - 1/2 < s le 3{{n - 1} over 2} + 3/2 - 1/2quad Rightarrow cr & Rightarrow quad leftlfloor {3{{n - 1} over 2} - 3/2 - 1/2} rightrfloor < s le leftlfloor {3{{n - 1} over 2} + 3/2 - 1/2} rightrfloor cr} $$ and in general, for a dimension $m$ $$ bbox[lightyellow] { eqalign{ & d - m < x_{,k} + y_{,j} + z_{,l} = s < dquad Rightarrow cr & Rightarrow quad m{{n - 1} over 2} - m/2 - 1/2 < s le m{{n - 1} over 2} + m/2 - 1/2quad Rightarrow cr & Rightarrow quad leftlfloor {{{mn - 1} over 2}} rightrfloor - m < s le leftlfloor {{{mn - 1} over 2}} rightrfloor cr} tag{3}}$$ which leads to $$ bbox[lightyellow] { N(n,m) = M_b left( {leftlfloor {{{mn - 1} over 2}} rightrfloor ,;n - 1,;m} right) - M_b left( {leftlfloor {{{mn - 1} over 2}} rightrfloor - m,;n - 1,;m} right) tag{4}}$$

The values for smaller $m$ and $n$ given by the formula are

which check against direct computation.

Finally, concerning the asymptotic for large $n$, we make the following considerations:

• from the sketch of the inequality above it is clear that is encompasses the $m$ most central bar of the $N_b$ histogram;
• for large values of $n$, being $N_b$ the convolution of three variables uniform over a wide stretch, it is evident that the central values will flatten out, and we can take for $N$ $m$ times the central value $$ N_b left( {leftlfloor {mleft( {n - 1} right)/2} rightrfloor ,n - 1,m} right) $$
• it is not easy to compute the peak value of $N_b$ in the general case (re. to this post ) but for $m=3$ it is quite straight: the number of integral points on the plane $x+y+z=s=3leftlfloor m(n-1)/2 rightrfloor $ correspond to those projected into the $x,y$ plane into the inequality $s-(n-1) le x+y le s-0$ as in this sketch

so that the maximum of $N_b$ equals the points in the central stripe as shown, for large $n$ (small unit squares) tending to the
continuous and thus giving $$ bbox[lightyellow] { N(n,3) approx {9 over 4}left( {n - 1} right)^2 tag{5}}$$

and in fact