We investigate sets of mutually orthogonal latin squares (MOLS) generated by cellular automata (CA) over finite fields. After introducing how a CA defined by a bipermutive local rule of diameter $d$ over an alphabet of $q$ elements generates a Latin square of order $q^{d-1}$, we study the conditions under which two CA generate a pair of orthogonal Latin squares. In particular, we prove that the Latin squares induced by two Linear Bipermutive CA (LBCA) over the finite field $\mathbb{F}_q$ are orthogonal if and only if the polynomials associated to their local rules are relatively prime. Next, we enumerate all such pairs of orthogonal Latin squares by counting the pairs of coprime monic polynomials with nonzero constant term and degree n over $\mathbb{F}_q$. Finally, we present a construction for families of MOLS based on LBCA, and prove that their cardinality corresponds to the maximum number of pairwise coprime polynomials with nonzero constant term. Although our construction does not yield all such families of MOLS, we show that the resulting lower bound is asymptotically close to their actual number.
Note: this publication is an extension of this exploratory paper presented at AUTOMATA 2016.