20 also determined solely by background cosmology. These background quantities should not bother us here. The scalar field EOM becomes In the MLAPM code the partial differential equation Eq. (31) is (and in our modified code Eq. (32) will also be) solved on discretized grid points, and as such we must develop the discretized versions of Eqs. (29 - 32) to be implemented into the code. But before going on to the discretization, we need to address a technical issue. As the potential is highly non- linear, in the high density regime the value of the scalar field √$_{κϕ}$ will be very close to 0, and this is potentially a disaster as during the numerical solution process the value of √$_{κϕ}$ might easily go into the forbidden region ϕ < 0 [25]. One way of solving this problem is to define χ = ¯ χe u in which ¯ χ is the background value of χ , as in [25]. Then the new variable u takes value in ( −∞ , ∞ ) so that e u is positive definite which ensures that χ > 0. However, since there are already exponentials of χ in the potential, this substitution will result terms involving exp [exp( u )], which could potentially magnify any numer- ical error in u . Instead, we can define a new variable u according to plus Eqs. (A2, A3). These equations will ultimately be used in the code. Among them, Eqs. (A2, A5) will use the value of u while Eq. (A3) solves for u . In order that these equations can be integrated into MLAPM, we need to discretize Eq. (A3) for the application of Newton-Gauss- Seidel iterations. in which we have used the fact that χ = log(1 + e $^{u}$) ⇒ ∇ χ = e u 1+ e u ∇ u , and moved all terms depending only on background cosmology (the source terms) to the right hand side. So, in terms of the new variable u , the set of equations used in the N -body code should be Then the Poisson equation becomes By this, u still takes value in ( −∞ , ∞ ), e u ∈ (0 , ∞ ) and thus e χ ∈ (1 , ∞ ) which ensures that χ is positive definite in numerical solutions. Besides, e βχ = [1 + e $^{u}$] β so that there will be no exponential-of-exponential terms, and the only exponential is what we have for the potential itself. β = − 1 above. To discretize Eq. (A3), let us define b ≡ e u 1+ e $_{u}$. The discretization involves writing down a discretion version of this equation on a uniform grid with grid spacing h . Suppose we require second order precision as is in the standard Poisson solver of MLAPM, then ∇ u in one di- mension can be written as ac 2 ( H$_{0}$B ) 2 ∇ · ( e u 1 + e u ∇ u ) = 3 γ Ω$_{CDM}$ ρ$_{c,}$$_{CDM}$ (1 + e $^{u}$) γ + 3 µβ Ω$_{V}$$_{0}$ a $^{3}$(1 + e $^{u}$) β [ 1 − (1 + e $^{u}$) β ] µ +1 − 3 γ Ω$_{CDM}$ e $_{γ}$√$_{κ}$$_{¯}$ ϕ − 3 µβ Ω$_{V}$$_{0}$ a $^{3}$e $_{β}$√$_{κ}$$_{¯}$ ϕ [ 1 − e $_{β}$√$_{κ}$$_{¯}$ ϕ ] m +1 (A3) ∇ $^{2}$Φ$_{c}$ = 3 2 Ω$_{CDM}$ [ ρ$_{c,}$$_{CDM}$ (1 + e $^{u}$) γ − e $_{γ}$√$_{κ}$$_{¯}$ ϕ ] + 3 2 Ω$_{B}$ ( ρ$_{c,}$$_{B}$ − 1) − 3Ω$_{V}$$_{0}$ a 3 [ 1 − (1 + e $^{u}$) β ] µ + 3 ¯ Ω$_{V}$ a $^{3}$, (A2) ∇ · ( b ∇ u ) → ( ∇ $^{h}$b$_{j}$ ) · ( ∇ $^{h}$u$_{j}$ ) + b$_{j}$ ∇ h $^{2}$u$_{j}$ = b$_{j}$$_{+1}$$_{/}$$_{2}$ − $^{b$_{j}$}$− 1 / 2 h u$_{j}$$_{+1}$ − $^{u$_{j}$}$− 1 2 h + b$_{j}$$_{+1}$$_{/}$$_{2}$ + $^{b$_{j}$}$− 1 / 2 2 u$_{j}$$_{+1}$ − 2 u$_{j}$ + $^{u$_{j}$}$− 1 h 2 = 1 h 2 [ b$_{j}$$_{+}$ 1 $_{2}$u$_{j}$$_{+1}$ − u$_{j}$ ( b$_{j}$$_{+}$ 1 2 + $^{b$_{j}$}$− 1 2 ) + $^{b$_{j}$}$− 1 $^{$_{2}$u$_{j}$}$− 1 ] (A7) ∇ u → ∇ $^{h}$u$_{j}$ = u$_{j}$$_{+1}$ − $^{u$_{j}$}$− 1 2 h (A6) d x c dt$_{c}$ = p c a $_{2}$, (A4) d p c dt$_{c}$ = − 1 a ∇ Φ$_{c}$ [ − c $^{2}$γ ( H$_{0}$B ) 2 e u 1 + e u ∇ u ] (A5) Appendix A: Discretized Equations for Our N -body Simulations e u + 1 = e $^{χ}$. (A1) [41] S. P. D. Gill, A. Knebe and B. K. Gibson, Mon. Not. R. Astron. Soc. 351 , 399 (2004). [42] S. Knollmann and A. Knebe, Astrophys. J. Suppl. 182 , 608 (2009). where we have defined ¯ Ω$_{V}$ = κV ( ¯ ϕ ) / 3 H 2 0 which is deter- mined by background cosmology, the quantity e $_{γ}$√$_{κ}$$_{¯}$ ϕ is where a subscript j means that the quantity is evaluated on the j -th point. Of course the generalization to three dimensions is straightforward. The factor b in ∇· ( b ∇ u ) makes this a standard variable coefficient problem. We need also discretize b , and do it in this way (again for one dimension):