We analyze ground state (GS) factorization in general arrays of spins si with XXZ couplings immersed in nonuniform fields. It is shown that an exceptionally degenerate set of completely separable symmetry-breaking GS's can arise for a wide range of field configurations, at a quantum critical point where all GS magnetization plateaus merge. Such configurations include alternating fields as well as zero bulk field solutions with edge fields only and intermediate solutions with zero field at specific sites, valid for d-dimensional arrays. The definite magnetization projected GS's at factorization can be analytically determined and depend only on the exchange anisotropies, exhibiting critical entanglement properties. We also show that some factorization compatible field configurations may result in field-induced frustration and nontrivial behavior at strong fields.One of the most remarkable phenomena arising in finite interacting spin systems is that of factorization. For particular values and orientations of the applied magnetic fields, the system possesses a completely separable exact ground state (GS) despite the strong couplings existing between the spins. The close relation between GS factorization and quantum phase transitions was first reported in [1] and has since been studied in various spin models [2][3][4][5][6][7][8][9][10][11][12], with general conditions for factorization discussed in [7] and [13]. Aside from some well known integrable cases [14-17], higher dimensional systems of arbitrary spin in general magnetic fields are not exactly solvable, so that exact factorization points and curves provide a useful insight into their GS structure.The XXZ model is an archetypal quantum spin system which has been widely studied to understand the properties of interacting many-body systems and their quantum phase transitions [18][19][20][21][22][23] Our aim here is to show that in finite XXZ systems of arbitrary spin under nonuniform fields, highly degenerate exactly separable symmetry-breaking GS's can arise for a wide range of field configurations in arrays of any dimension, at an outstanding critical point where all magnetization plateaus merge and entanglement reaches full range. The Pokrovsky-Talapov (PT)-type transition in a spin-1/2 chain in an alternating field [20] is shown to correspond to this factorization. Magnetization phase diagrams, showing non trivial behavior at strong fields, and pair entanglement profiles for distinct factorization compatible field configurations are presented, together with analytic results for definite magnetization GS's.We consider an array of N spins s i interacting through XXZ couplings and immersed in a general nonuniform magnetic field along the z axis. The Hamiltonian readswith h i , S µ i the field and spin components at site i and J ij , J ij z the exchange coupling strengths. Since H commutes with the total spin component S z = i S z i , its eigenstates can be characterized by their total magnetization M along z. The exact GS will then exhibit definite M plateaus a...