Calculations based on Chapter 5: Radial Distributions

In this notebook, we will attempt to solve the mystery raised in the previous notebook: Why do the core-ionization energies for an isonuclear series exhibit an approximately piecewise linear dependence on the number of electrons? This is surprising, because the effect on the wave function of adding electrons is complicated. How can it lead to such a simple result?

In order to understand this piecewise linear dependence, we will calculate the potential experienced by an electron in the 1s-shell due to the other electrons. This potential is a function of the number of electrons, N. We hold Z constant and increase N. To calculate the potential, we first calculate the radial electron density distribution as a function of N.

Initialization

Get["Sturmian`"]; (* We include this for hwf[], the definition of
                      the one-electron hydrogen-like wave functions. *)
 Get["Sturmian`binary`"];  (* For fast repulsionMatrix[] *)
 Get["Sturmian`M`pauli`"];
 Off[General::spell];

Radial distribution of electron density

In Section 2 of Chapter 5, we show analytically, that if the electrons in the n==2 shell were confined to a thin shell, then the potential experienced by an electron inside that shell would be proportional to N for 2≤N≤10. This crude model explains qualitatively the lineariarity of the core ionization energies for this range of N. However, we would like to explain the linearity using a realistic, calculated charge distribution.

We begin by noting, that the electron density of a multiconfigurational wave function Ψ==(∑_ν)(Subscript[Φ, ν])(Subscript[B, ν]) is given by

({{ρ(x), ==, N∫ds∫d(Subscript[τ, 2])⋯∫d(Subscript[τ, N])(SuperStar[Ψ])(x,(Subscript[x, 2]),…,(Subscript[x, N]))Ψ(x,(Subscript[x, 2]),…,(Subscript[x, N]))}, {, ==, (∑_(ν',ν))(Subscript[B, ν']^( * ))(Subscript[ρ, ν'ν])(x)(Subscript[B, ν])    (x∈(R^3))}})

where

(Subscript[ρ, ν'ν])(x)==N∫ds∫d(Subscript[τ, 2])⋯∫d(Subscript[τ, N])(Subscript[Φ, ν']^( * ))(x,(Subscript[x, 2]),…,(Subscript[x, N]))(Subscript[Φ, ν])(x,(Subscript[x, 2]),…,(Subscript[x, N]))

Here dτ indicates integration over space and spin coordinates, while ds indicates spin alone. In the case when (Subscript[ℛ, ν])==(Subscript[ℛ, ν']), the one-electron configurations (Subscript[χ, μ]) are orthonormal, and we get

({{(Subscript[ρ, ν'ν])((Subscript[x, 1])), ==, ({{0,  for , |ν'\ν|≥2}, {, , }, {(Subscript[χ, μ']^( * ))((Subscript[x, 1]))(Subscript[χ, μ])((Subscript[x, 1])),  for , ν'\ν=={μ} and ν\ν'=={μ'}}, {, , }, {(∑_(i==1)^N)|(Subscript[χ, (Subscript[μ, i])])((Subscript[x, 1]))((|)^2),  for , ν'==ν}})}})

such that the full density becomes

ρ(x)==(∑_ν)|(Subscript[B, ν])((|)^2)(∑_(j==1)^N)|(Subscript[χ, (Subscript[μ, j])])(x)((|)^2)+(∑_(|ν'\ν|==1))(Subscript[χ, μ'])(SuperStar[(x)])(Subscript[χ, μ])(x)

To find the spherically averaged distribution function ρ(r), we must then average over solid angle; but because of the orthonormality relations for spherical harmonics, this simply yields

ρ(r)==(1/(4π))∫ρ(x)dΩ==(1/(4π))(∑_ν)|(Subscript[B, ν])((|)^2)(∑_(j==1)^N)|(Subscript[R, (Subscript[μ, j])])(r)((|)^2)+(1/(4π))(∑_(|ν'\ν|==1))(Subscript[R, μ'])(r)(Subscript[R, μ])(r)

The radial distribution function, which will be shown in the plots, is (r^2)ρ(r). This function shows the shell structure more clearly than the spherically averaged density distribution, which is sharply peaked at r==0 and only decays with r.

The parameter p below is a predicate on the μ’s, deciding for each one whether or not to include it in the sum. This allows us to look at the contributions to the density from particular shells. If no p is provided, then everything is included.

(* The spherically averaged density distribution (Subscript[ρ, ν'ν])(r) *)
partialdensity[νp_, ν_, Q_, p_] := 
  Module[{q, ns, ls, Nel = Length[ν], n, l, m, ms, np, lp, mp, 
    msp,diff,contributions},
  
   contributions = Select[ν,p];
   ns = #[[1]] & /@ contributions;
   ls = #[[2]] & /@ contributions;

   diff = Complement[νp, ν];
   (4 Pi)^(-1)*
    Switch[Length[diff],
     (* Diagonal elements *)
     0, Sum[(hwf[ns[[i]], ls[[i]], Q, r]*Exp[-(Q/ns[[i]])*r])^2, {i, 1, Length[ns]}],

     (* ν and νp differ by a single μ *)
     1, ({n, l, m, ms} = diff[[1]]; {np, lp, mp, msp} = Complement[ν, νp][[1]];
         If[p[{n,l,m,ms}] && p[{np,lp,mp,msp}],
            hwf[n, l, Q, r]*hwf[np, lp, Q, r]*Exp[-(Q/n + Q/np)*r], 0]),

     (* νp and ν differ by more than a single μ *)
     _,0
     ]
   ];

partialdensity[νp_, ν_, Q_] := partialdensity[νp, ν, Q, True&];

(* Spherically averaged density distribution of a multiconfigurational state 
   within which (Subscript[ℛ, ν]) is constant *)
fulldensity[B_, Q_, basis_, p_] := 
 Sum[B[[i]]*B[[j]]*partialdensity[basis[[i]], basis[[j]], Q, p], 
     {i, Length[basis]}, {j, Length[basis]}]

fulldensity[B_, Q_, basis_] := fulldensity[B, Q, basis, True&];

We shall start by recomputing the bases and inter-electron repulsion matrices that we calculated in the previous notebook.

{n1, n2, n3} = Table[oneelectronstates[n, n], {n, 3}];
n3sp = Select[n3, #[[2]] <= 1 &];
n2shells = Table[Join[n1, #] & /@ binomialchoices[n2, N2], {N2, 0, 8}];
n3shells = Table[Join[Join[n1, n2], #] & /@ binomialchoices[n3sp, N3], {N3, 1, 8}];

rblocks0 = Join[{{{{1, 0, 0, -1}}, {{1, 0, 0, 1}}}}, n2shells, n3shells];
rblocksi = Drop[#, 1] & /@ # & /@ rblocks0;  

Ts0 = repulsionMatrix /@ rblocks0;
Tsi = repulsionMatrix /@ rblocksi;

Then we will recompute the eigenvalues 
 and eigenvectors cs when 1≤N≤18:

eigen0 = Reverse /@ # & /@ Eigensystem /@ Ts0;
  λ0 = #[[1, 1]] & /@ eigen0;
  cs0 = Chop[#[[2, 1]]] & /@ eigen0;

  eigeni = Reverse /@ # & /@ Eigensystem /@ Tsi;
  λi = #[[1, 1]] & /@ eigeni;
  csi = Chop[#[[2, 1]]] & /@ eigeni;

In generalized Sturmian calculations, each configuration is characterized by an effective charge Q==((Subscript[p, κ])/(Subscript[ℛ, ν])). Since in this calculation we are working with the large-Z approximation (i.e. within a fixed (Subscript[ℛ, ν])-block), we have (Subscript[p, κ])==(Subscript[Zℛ, ν])+(Subscript[λ, k]) and so

Q==Z+((Subscript[λ, κ])/(Subscript[ℛ, ν]))

Thus within the (Subscript[ℛ, ν])-block, all configurations contributing to a particular state (Subscript[Ψ, κ]) will have the same effective charge Q. We also recall that

((Subscript[ℛ, ν])^2)==((Subscript[N, 1])/(1^2))+((Subscript[N, 2])/(2^2))+((Subscript[N, 3])/(3^2))+⋯

where (Subscript[N, j]) is the number of μ∈ν that have n==j. Thus, when 1≤N≤18, the (Subscript[ℛ, ν]) of the “ground state” ℛ-block and of the “core ionized ground state” ℛ-block will be

Rν0Sqr[Nel_] := Module[
   {N1 = Min[Nel, 2], N2 = Max[0, Min[Nel - 2, 8]], N3 = Max[0, Min[Nel - 10, 28]]},
   N1 + N2/4 + N3/9
 ];
 RνiSqr[Nel_] := Rν0Sqr[Nel] - 1;

and the corresponding Q’s are:

Qs0 = Table[Z + λ0[[Nel]]/ Sqrt[Rν0Sqr[Nel]], {Nel, 1, 18}];
 Qsi = Join[{0},Table[Z + λi[[Nel]]/ Sqrt[RνiSqr[Nel]], {Nel, 2, 18}]];

To get a little intuition, let us calculate the radial distribution of eletron density for a few states. The following computes and plots (r^2)ρ(r) for the large-Z approximation to the ground state of oxygen (N==8, Z==8).

{Nel,Z0} = {8,8}; (* Number of electrons, Nuclear charge *)
wf = cs0[[Nel]] (* Coefficients for ground state *)

basis = rblocks0[[Nel]];
ρ = fulldensity[wf, Qs0[[Nel]] /. Z -> Z0, basis] // Simplify;
Ogroundplot = Plot[Evaluate[r^2 ρ], {r, 0, 2}, PlotRange -> All, 
                   AxesLabel -> {r,"$r^2$ ρ(r) = " r^2 ρ}, 
                   ImageSize -> 700]

Next, we compute the radial distribution function for the core-ionized ground state of oxygen. The reader is encouraged to inspect the density distributions for other atoms and ions (remembering that the computed roots restrict N to 2,…,10).

{Nel,Z0} = {8,8};   (* Number of electrons, Nuclear charge *)
wf = csi[[Nel]] (* Coefficients for core-ionized ground state *)

basis = rblocksi[[Nel]];
ρ = fulldensity[wf, Qsi[[Nel]] /. Z -> Z0, basis] // Simplify;
Oplot = Show[
    Ogroundplot,
    Plot[Evaluate[r^2 ρ], {r, 0, 2}, PlotRange -> All, 
                   AxesLabel -> {r,r^2 "ρ(r) = " r^2 ρ}, 
                   ImageSize -> 700,PlotStyle->Red]
]

It can be seen in the figure, that the contribution from the 1s-shell is reduced by half when one of the electrons is removed from it; and the n==2 is drawn in slightly towards the nucleus, because the nucleus is less screened when an electron is removed from the inner shell.

In the previous notebook, we encountered the remarkable fact that, for isonuclear series, the core-ionization energy obeys an approximately piecewise linear relationship. To understand this we will now generate first the electronic densities and then the electronic potentials for an isonuclear series; that is to say, we hold Z constant and generate the densities and potentials corresponding to various values of N. When a core electron is lifted from its orbital to an ionized state, it experiences a potential that is partly due to the attraction of the nucleus and partly due to the repulsion of all the other electrons. For an isonuclear series, the first term will be constant, while the second term, which we will calculate below, will vary with the number of electrons N.

Radial Density Distributions for an Isonuclear Series

In this section, we generate the electronic densities for an isonuclear series and from these, in the next section, we will calculate the electronic potentials. The following computes the contribution to the density from the n==2 shell (dn2shell) and the full density (dallshells) for the core-ionized ground states when the nuclear charge is fixed to Z==18 while N assumes the range from 2 to 10.

Z0 = 18;
dn2shell = Table[fulldensity[Chop[csi[[Nel]]], Qsi[[Nel]] /. Z -> Z0, 
                  rblocksi[[Nel]],#[[1]]==2&], {Nel, 2, 10}] // Simplify;

dallshells = Table[fulldensity[Chop[csi[[Nel]]], Qsi[[Nel]] /. Z -> Z0, 
                    rblocksi[[Nel]]], {Nel, 2, 10}] // Simplify;

We can take a look at the contributions from the n==2 shell:

Table[{"N = ", n, " ρ(r) = ", dn2shell[[n - 1]]}, {n, 2, 10}] // TableForm

The following visualizes first the contributions to the density from the n==2 shell, then the full densities for the core-ionized states:

Plot[Evaluate[r^2*dn2shell], {r, 0, 2}, PlotRange -> All]
Plot[Evaluate[r^2*dallshells], {r, 0, 2}, PlotRange -> All]

For comparison, densities for the ground states of the isonuclear series:

dallshells0 = Table[fulldensity[Chop[cs0[[Nel]]], Qs0[[Nel]] /. Z -> Z0,  
                     rblocks0[[Nel]]], {Nel, 2, 10}] // Simplify;

Plot[Evaluate[r^2*dallshells0], {r, 0, 2}, PlotRange -> All]

We can check that the radial density function integrates to the number of electrons that contribute to it:

Integrate[4 Pi r^2 * dallshells, {r, 0, Infinity}]

Integrate[4 Pi r^2 * dallshells0, {r, 0, Infinity}]

Potential experienced by ionized core-electron

The electronic potential φ(x), i.e. the contribution to the potential due to the electrons, is related to the electronic density ρ(x) by

φ(x)==(∫_(R^3))dx'((ρ(x'))/(x-x'))

Starting instead with the electronic radial density function of a given state, we can thus construct the spherically averaged electronic potential at a given radius as:

φ(r)==4π(∫_0^r)dr'(r^(′2))ρ(r')+4π(∫_r^∞)dr'r'ρ(r')

potential[c_, Q_, basis_, r0_, p_] := Module[{ρ, int1, int2},
      ρ = fulldensity[c, Q, basis, p] /. r -> rp // Simplify;
      
   int1 = Integrate[rp^2* ρ, {rp, 0, r0}];
   int2 = Integrate[rp ρ, {rp, r0, Infinity}];
       
   If[r0 > 0, (4 Pi/r0)*int1, 0] + 4 Pi*int2
];
potential[c_, Q_, basis_, r0_] := potential[c, Q, basis, r0, True &];

In this section, using our calculated core-ionized ground states, we will compute and visualize the contribution to the potential by the electrons in the n==2 shell, when the nuclear charge is fixed and 2≤N≤10.

Z0 = 10;
potsi = Table[potential[csi[[Nel]], Qsi[[Nel]] /. Z -> Z0, 
               rblocksi[[Nel]],r,#[[1]]==2&], {Nel, 2, 10}];

As a sanity check on the calculated potentials, let us inspect them visually:

potplots  = Plot[Evaluate[potsi], {r, 0, 2}, AxesLabel -> {r}, PlotRange->All]

We observe that the spherically averaged electronic potentials fall in two regions: an inner region, where the beahviour is non-trivial, and an outer region in which the potentials decay monotonically. This is because, outside of an atom or atomic ion, an electron experiences its electronic potential approximately as if the electronic charges were all concentrated at its center.Of course, the atomic radius is not well-defined: Electrons may move anywhere in space. However, for bound states, the overwhelming bulk of the probability is concentrated within a bounded region, which we may call the “atomic region”.

The following shows that the electronic potentials, outside of the atomic region, become approximately proportional to (N/r). The vertical line in the second graph shows the atomic radius of Ne: 38pm ≃ 0.72 Bohr.

Neline = Graphics[Line[{{.72, 0}, {.72, 10}}]];
rpotplots  =Plot[Evaluate[r*potsi], {r, 0, 2}, AxesLabel -> {r}]
rpotplotsN = Show[
  	Plot[Evaluate[Table[r/N2*potsi[[N2 + 1]], {N2, 1, 8}]], {r, 0, 3}, 
    	     AxesLabel -> {r, r/Subscript[N, 2] φ[r]}, 
     	     PlotRange -> {{0, 3}, {0, 1}}],
        Neline
]

Having calculated the potentials, we can make a plot of the potential at the nucleus as a function of the number of electrons N for an isonuclear series. Here we consider the isonuclear series with Z==18. The table 
0 generates the potential from the multiconfigurational wave function obtained by diagonalizing T'. The reader can find corresponding results for other isonuclear series by changing Z0.

Z0 = 18;
φ0 = Table[potential[cs0[[Nel]], Qs0[[Nel]] /. Z -> Z0, rblocks0[[Nel]], 0], {Nel, 1, 18}]

ListPlot[φ0,AxesOrigin->{0,0},AxesLabel -> {"N","φ(0)"}]

To see the effect of interelectron repulsion at the nucleus, we generate this potential as it would be with T'==0. In this case, we do not need a multiconfigurational wave function; we used the single configuration chosen according to Hund’s rules. The code below implements the special case of computing the density and electronic potential from the lowest-energy single-configuration wave function.

lowconfiguration = {{1, 0, 0, 1}, {1, 0, 0, -1}, {2, 0, 0, 1}, {2, 0, 0, -1}, 
  {2, 1, -1, 1}, {2, 1, 0, 1}, {2, 1, 1, 1}, {2, 1, -1, -1}, {2, 1, 0, -1}, 
  {2, 1, 1, -1}, {3, 0, 0, 1}, {3, 0, 0, -1}, {3, 1, -1, 1}, {3, 1, 0, 1}, 
  {3, 1, 1, 1}, {3,  1, -1, -1}, {3, 1, 0, -1}, {3, 1, 1, -1}, {3, 2, -2, 1}, 
  {3, 2, -1, 1}, {3, 2, 0, 1}, {3, 2, 1, 1}, {3, 2, 2, 1}, {3, 2, -2, -1}, 
  {3, 2, -1, -1}, {3, 2, 0, -1}, {3, 2, 1, -1}, {3, 2, 2, -1}};

lowconfig[n_] := Take[lowconfiguration, n];

(* Radial density distribution of lowest energy N-electron configuration *)
lowdensity[Nel_, Q_] := partialdensity[lowconfig[Nel], lowconfig[Nel], Q];

lowpotential[n_, Q_, r0_] := Module[{ρ},
  ρ = Simplify[lowdensity[n, Q] /. r -> rp];
  

   If[r0 > 0, (4 Pi/r0)*Integrate[rp^2 ρ, {rp, 0, r0}], 0]
    + 4 Pi*Integrate[rp ρ, {rp, r0, Infinity}]
  ]

The table φ00 neglects interelectron repulsion as described in the previous paragraph. We plot the computed values together with φ0.

φ00  = Table[lowpotential[Nel, Z0, 0], {Nel, 18}]

piecewisedots = Show[
  ListPlot[φ00, PlotStyle -> {{RGBColor[0, 0, .4], PointSize[.007]}}],
  ListPlot[φ0, PlotStyle -> {{RGBColor[.4, 0, 0], PointSize[.007]}}],
  
  PlotRange -> {{0, 21.5}, {0, 105}}, AxesOrigin -> {0, 0}, 
  AxesLabel -> {"N", "φ(0)"},
  ImageSize -> 500
  ]

In this figure, the upper dots correspond to densities generated neglecting interelectron repulsion. In the lower set of dots, interelectron repulsion is included. In the upper set, the points corresponding to the filling of the n==2 and n==3 shells, respectively, can be fitted exactly by straight lines, as shown below

Fit[Take[φ00, {3, 10}], {1, (Subscript[N, 2]), (Subscript[N, 2])^2, (Subscript[N, 2])^3}, (Subscript[N, 2])] // Chop
Fit[Take[φ00, {11, 18}], {1, (Subscript[N, 3]),(Subscript[N, 3])^2, (Subscript[N, 3])^3}, (Subscript[N, 3])] // Chop

The fact that the contribution from each shell is exactly proportional to the number of electrons in the shell when interelectron repulsion is neglected, suggests to us that there is an exact relationship that could be discovered, and that this relationship will also inform us as to the approximate proportionality when interelectron repulsion is included.

Let us now investigate this question. Within the framework of the large-Z approximation we have

∫dτ (Subscript[Ψ, κ]^( * ))(x)(Subscript[V, 0])(x)(Subscript[Ψ, κ])(x)==(∑_(ν'))(∑_ν)(Subscript[B, ν'κ]^( * ))(Subscript[B, νκ])∫dτ (Subscript[Φ, ν']^( * ))(x)(Subscript[V, 0])(x)(Subscript[Φ, ν])(x)==-((Subscript[p, κ]^2)/(Subscript[β, ν]))(∑_ν)|(Subscript[B, νκ])((|)^2)

In the second step above, we make use of the potential weighted orthonormality relation. Further, since (∑_ν)|(Subscript[B, νκ])((|)^2)==1, this equation reduces to

∫dτ (Subscript[Ψ, κ]^( * ))(x)(Subscript[V, 0])(x)(Subscript[Ψ, κ])(x)==-((Subscript[p, κ]^2)/(Subscript[β, ν]))==-(Subscript[p, κ])(Subscript[Zℛ, ν])

This result can be used to express the electronic potential at the nucleus in a very simple form. Combining the definitions of φ(x) and ρ(x), we obtain

φ(0)==N∫dτ(1/(Subscript[x, 1]))(Subscript[Ψ, κ]^( * ))(x)(Subscript[Ψ, κ])(x)

From the definition of (Subscript[V, 0]) and from the fact that each term in the sum (Subscript[V, 0])==-(∑_(i==1)^N)(Z/(Subscript[r, i])) gives the same contribution, we have

φ(0)==-(1/Z)∫d(Subscript[τΨ, κ]^( * ))(x)(Subscript[V, 0])(x)(Subscript[Ψ, κ])(x)

Thus we obtain the extremely simple result:

φ(0)==(Subscript[p, κ])(Subscript[ℛ, ν])

which can alternatively be written in the form:

φ(0)==((Subscript[Zℛ, ν])^2)-|(Subscript[λ, κ])|(Subscript[ℛ, ν])

or in a third form:

φ(0)==(Subscript[Q, ν])((Subscript[ℛ, ν])^2)

since (Subscript[Q, ν])==Z-|(Subscript[λ, κ])|/(Subscript[ℛ, ν]). From this, it follows that for an isonuclear series, the electronic potential at the nucleus depends on N in an approximately piecewise linear way. For example, let us consider the isonuclear series where Z==18. Keeping the nuclear charge Z constant at this value, we begin to add electrons. For the ground state we have:

((Subscript[ℛ, ν])^2)≡(1/(Subscript[n, 1]^2))+(1/(Subscript[n, 2]^2))+⋯+(1/(Subscript[n, N]^2))==({{(N/1),     N≤2}, {, }, {(2/1)+((N-2)/4),     2≤N≤10}, {, }, {(2/1)+(8/4)+((N-10)/9),     10≤N≤18}})

We can plot the closed form expressions simultaneously with the electronic potentials at the nucleus, laboriously derived from calculated densities.

φ0closed   = Table[Qs0[[Nel]]*Rν0Sqr[Nel] /. Z -> Z0, {Nel, 1, 18}]

labels = Graphics /@ {
  Text["ZN", {4.2, Z0*4.2}], Text["2Z+((N-2)/4)",  {13.3, 2 Z0 + Z0 (13.8 - 2)/4}],
  Text["4Z+((N-10)/9)", {20, 4 Z0 + Z0 (20 - 10)/9 + 2}], Text["(Subscript[Q, ν])((Subscript[ℛ, ν])^2)", {19, φ0[[18]] + 2}]
};

piecewise = Show[
  ListPlot[φ0closed, PlotStyle -> {{Red}}, Joined -> True, 
   PlotRange -> {{0, 21.5}, {0, 105}}, AxesOrigin -> {0, 0}, 
   AxesLabel -> {"N", "φ(0)"}],
  
  labels,
  ListPlot[φ0, PlotStyle -> {{RGBColor[.4, 0, 0], PointSize[.007]}}],
  
  Plot[Z0 Nel, {Nel, 0, 4}, PlotStyle -> Dashed],
  Plot[2 Z0 + Z0 (Nel - 2)/2^2, {Nel, 1, 13}, PlotStyle -> Dashed],
  Plot[4 Z0 + Z0 (Nel - 10)/3^2, {Nel, 8, 19.5}, PlotStyle -> Dashed],
  ListPlot[φ00, PlotStyle -> {{RGBColor[0, 0, .4], PointSize[.007]}}],
  
  ImageSize -> 500
  ]

Qualitative Connection to Core-Ionization Energies

If the effect of the electrons in the shells with n≥2 were neglected, the core-ionization energy would be that of helium. A rough estimate of the effect of the outer shells can be obtained by subtracting the potential at the nucleus due to them from the ionization energy of helium. That is, the core-ionization will be the (N-independent) ionization energy of helium minus φ(0) due to the n==2 and higher shells. This approximation corresponds to the assumption that the densities of the outer shells do not overlap the density of the 1s-shell.

In reality, the density distribution of the 1s-electron is localized near to the nucleus, but overlaps the densities from the shells of higher n. This can be seen in the following figures, where the contributions to the densities from the n==1, n==2 and n==3 shells are shown seperately - alongside the full density distributions - for the Z==18 isonuclear series.

{Z0, Nmax} = {18, 18};

colours = Table[RGBColor[1 - N[i/Nmax], 1/i, N[i/(2 Nmax)]], {i, Nmax }];
fills = Table[i -> {{i + 1}, colours[[i]]}, {i, Nmax}];
range = {0, 2.1};
style = Opacity[.4];

dshell[n_] := Join[{0},Table[fulldensity[Chop[cs0[[Nel]]], Qs0[[Nel]] /. Z -> Z0, 
                    rblocks0[[Nel]], #[[1]] == n &], {Nel, 1, Nmax}] ];
shells = dshell /@ {1, 2, 3};

n1shell   = Plot[Evaluate[r^2*shells[[1]]], {r, 0, 2}, PlotRange -> range,Filling -> fills, PlotStyle -> style]
n2shell   = Plot[Evaluate[r^2*shells[[2]]], {r, 0, 2}, PlotRange -> range,Filling -> fills, PlotStyle -> style]
n3shell   = Plot[Evaluate[r^2*shells[[3]]], {r, 0, 2}, PlotRange -> range,Filling -> fills, PlotStyle -> style]
allshells = Plot[Evaluate[r^2*Total[shells]], {r, 0, 2}, PlotRange -> range,Filling -> fills, PlotStyle -> style]

Although the approximation of assuming the n==1 shell to be disjoint from shells with higher n does not completely hold, the piecewise-linear dependence of φ(0) on the number N of electrons for isonuclear series can give us a qualitative understanding of the piecewise linear dependence on N of the core-ionization energies.