weird assumption of one-way random effects ANOVA $$Y_{ij}  =  μ \, +\, A_i \, +\,ε_{ij}$$
In Ch5 (random effects one-way ANOVA) of my textbook, it mentioned that $A_i$ (see the model above) is assumed to be $\sim N( 0, \text{constant_variance} )$
However, I can't find on the internet anything about this assumption. For example, look at the link: http://influentialpoints.com/Training/one-way_random_effects_anova-principles-properties-assumptions.htm 
My question is:


*

*Does this assumption have to hold to make the valid statistical inference?

*continuation of the first question, if yes, how can we check this assumption? I know we can check the assumptions of εij by looking into residuals(e.g. homogeneity, normal, etc.. ), but I don't know how to check the $A_i \sim N( 0, \text{constant_variance} )$ assumption.
p.s. this is my textbook : https://www.amazon.com/Design-Experiments-Statistical-Principles-Research/dp/0534368344

update:
I do some math and it seems that the assumption must hold !!!
Here is why:
When calculating the F-statistic(MSTR/MSE), the Expectation of Ai must be zero or otherwise we won't get correct MSTR(SSTR/df). please see the link below for calculating SSTR 
https://drive.google.com/file/d/1m2iMOU4CD-TUOdinNe3t-ylNdFMlzf2r/view?usp=sharing
But.. how can we check this assumption ??
 A: I am coming at this from a mixed models regression framework, not ANOVA necessarily, but it still applies (if it doesn't apply, I'm sure someone will let me know). But my terminology may be a little different.
You can think of the $\mu$ value "centering" the distribution of $a_{i}$.
Let's imagine you and I are in a group of friends that plays an online video game together, like Overwatch. Players accrue points throughout a match. Imagine I were to record how many points each of us scored over 30 games, and there were 10 people in this group of friends. If we average across all recorded games, let's say that the average points per match was 1500. That means $\mu = 1500$.
Now, each person has their own average, too. Let's say that the 10 personal averages differ normally from the $\mu$ with a standard deviation of 200. That would mean $a_i \sim N(0, \sigma_a^2)$ where $\sigma = 200$. Note that the mean here has to be zero, because it is how much each person differs from the grand mean. If it wasn't zero, then the grand mean would be different.
Lastly, each of us 10 people would have 30 measurements each. So each measurement $j$ can differ from the personal mean that is $\mu + a_i$. We assume these follow a normal distribution with some variance: $\epsilon_i \sim N(0, \sigma_a^2)$. For instance, that standard deviation could be 100.
Put crudely, the $a$ part is how far a person is away from the overall mean, whereas the $\epsilon$ part is how far each observation is away from the person's mean. Since both of these represent differences from a mean, they, by definition, must have a mean of zero.
How would this look? And how can we test for normality? See some R code below. This simulates a "one-way random-effects model," which is basically just an intercept-only multilevel model with random intercepts.
library(lme4)
library(dplyr)
# set parameters ---------------------------------------------------------------
n_i <- 10           # number of people
n_obs <- 30         # number of observations per person
n_j <- n_i * n_obs  # number of total observations
mu <- 1500          # defining the grand mean

# simulate ---------------------------------------------------------------------
set.seed(1839) # set seed for reproducibility

# for each person, simulate how far their mean is from grand mean:
a_i <- rnorm(n_i, 0, 200)

# simulate data set:
dat <- tibble(
  mu = mu,                       # everyone comes from same grand mean
  id = rep(1:n_i, each = n_obs), # assign ids to observations
  a = a_i[id],                   # assign the personal deviation
  e = rnorm(n_j, 0, 100),        # assign observation-level deviation
  y = mu + a + e                 # simulate dependent variable
)

# fit model --------------------------------------------------------------------
mod <- lmer(y ~ 1 + (1 | id), dat)
summary(mod)

This summary returns:
Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 33262    182.4   
 Residual              9525     97.6   
Number of obs: 300, groups:  id, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  1436.11      57.95   24.78

A few things to note that the simulation worked: The overall estimate is 1436, not too far off from the 1500 we set in the code above. The standard deviation around this, by id (or $a$), is 182, not far from 200 we set in the code above. The standard deviation of residual is 98, close to the 100 we set in the code above.
# check assumptions ------------------------------------------------------------
# are residuals normally distributed?
qqnorm(residuals(mod)); qqline(residuals(mod))

I'm checking assumptions using a QQ-plot. The last line above returns:

Which is the plot for $\epsilon$. We can also look at $a$, which is more difficult to tell, since there are only 10 people in this fictitious example:
# are a normally distributed around zero?
qqnorm(ranef(mod)$id[[1]]); qqline(ranef(mod)$id[[1]])


In short, you can extract the $a$ parts—also known as the random intercepts in the R package I'm using—and see if they are normally-distributed, just like you would with residuals.
