I am physicist from the first  education, so it is apparently my field. As it follows from my experience in 3D FEM testing with application to the magnetic field calculation there is a problem with equation $\nabla \times (\nu \nabla \times \vec {A})=\vec {j}$. Therefore we prefer another form of this equation, for example, this one $\nabla \nu \times (\nabla \times \vec {A})+\nu \nabla \times \nabla \times\vec {A} =\vec {j}$ (ungauged form).
Then if we have Coulomb gauge $\nabla.\vec {A}$, it automatically turns to $\nabla \nu \times (\nabla \times \vec {A})-\nu \nabla ^2\vec {A} =\vec {j}$ (Coulomb gauge). Now we can compare two form using `mesh` from Tim Laska answer (thanks to him), and function `appro` from  xzczd answer (thanks to him as well). Let check Coulomb gauge first:

    u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro = 
     With[{k = 1. 10^4}, Tanh[k #]/2 + 1/2 &];
    \[Nu]1 = Simplify`PWToUnitStep@
        PiecewiseExpand@If[x <= a && y <= a && z <= a, 1/(\[Mu]r), 1] /. 
       UnitStep -> 
        appro;(*permeability depending on iron cube in mesh*)\
    \[CapitalGamma]d = {DirichletCondition[ux[x, y, z] == 0, y == 0], 
       DirichletCondition[ux[x, y, z] == -fluxDensity*b/2, y == b], 
       DirichletCondition[uy[x, y, z] == 0, x == 0], 
       DirichletCondition[uy[x, y, z] == fluxDensity*b/2, x == b], 
       DirichletCondition[uz[x, y, z] == 0, 
        y == b || y == 0 || x == 0 || x == b || z == 0 || z == b]};
    \[CapitalGamma]n = {0, 0, 0};
    op7 = Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] - \[Nu]1*
       Laplacian[u, {x, y, z}];(*Coulomb gauged*){mvpAx, mvpAy, mvpAz} = 
     NDSolveValue[{op7 == {0, 0, 0}, \[CapitalGamma]d}, {ux, uy, 
       uz}, {x, y, z} \[Element] mesh]; 
Visualization
[![Figure 1][1]][1]
 
Now let check ungauged form

    u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro = 
     With[{k = 2. 10^4}, ArcTan[k #]/Pi + 1/2 &];
    \[Nu]1 = Simplify`PWToUnitStep@
        PiecewiseExpand@If[x <= a && y <= a && z <= a, 1/(\[Mu]r), 1] /. 
       UnitStep -> 
        appro;(*permeability depending on iron cube in mesh*)\
    \[CapitalGamma]d = {DirichletCondition[ux[x, y, z] == 0, y == 0], 
       DirichletCondition[ux[x, y, z] == -fluxDensity*b/2, y == b], 
       DirichletCondition[uy[x, y, z] == 0, x == 0], 
       DirichletCondition[uy[x, y, z] == fluxDensity*b/2, x == b], 
       DirichletCondition[uz[x, y, z] == 0, 
        y == b || y == 0 || x == 0 || x == b || z == 0 || z == b]};
    \[CapitalGamma]n = {0, 0, 0};
    op7 = Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] - \[Nu]1*
       Laplacian[u, {x, y, z}]; op8 = 
     Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] + \[Nu]1*
       Curl[Curl[u, {x, y, z}], {x, y, z}];(*Coulomb gauged*){mvpAx, 
      mvpAy, mvpAz} = 
     NDSolveValue[{op8 == {0, 0, 0}, \[CapitalGamma]d}, {ux, uy, 
       uz}, {x, y, z} \[Element] mesh]; 
[![Figure 2][2]][2]
   
It looks reasonable but pay attention how we play with `k` and with `Tanh[]`(Coulomb gauge) and `ArcTan[]` (ungauged form). 

  [1]: https://i.sstatic.net/0AbkX.png
  [2]: https://i.sstatic.net/lkfT4.png