Contents 1 From Atoms to Solids 1.1 Packing Ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Madelung Energy of 2D Bipartite Lattices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Simultaneous Energy and Momentum Eigenstates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Bond Formation in 1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Bonding in Water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Excitations of the Free-Electron Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Attractive Component of Lennard-Jones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Approximate Morse Potential Eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 Kronig-Penney Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5 5 12 14 15 20 20 20 22 26 27
2 Electrons in Crystals: Translational Periodicity 35 2.1 Reciprocal Lattice Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.2 Periodic Part of the Single-Particle Wavefunction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.3 Effective Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.4 Fermi Sphere of a 2D Square Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.5 Band Structures of FCC Crystals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.6 Orthogonality of Bloch States . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.7 Band Width and Number of Nearest Neighbors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.8 Bloch Theorem for Many-Body Wavefunctions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 2.9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 2.10 1D and 2D free-electron DOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 2.11 Critical Points in 1D and 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.12 DOS at Dirac point in 1D, 2D, 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.13 Graphene Band-Structure in Tight-Binding Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.14 Tight-Binding Model of Copper(II) Peroxide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3 Symmetries Beyond Translational Periodicity 71 3.1 Proper Rotations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.2 Symmetries of the 2D Square Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.3 Symmetries of the 2D Honeycomb Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.4 Special k-Points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.5 Point-Group Symmetries of the Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.6 Optical Transitions of NV Center in Diamond . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4 From Many Particles to the Single-Particle Picture 89 4.1 Hydrogen Molecule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.2 Many-Body States of the NV Center in Diamond . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.3 Hartree-Fock Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4.4 Single Particle Spectrum of Hartree-Fock Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.5 Bulk Modulus of Free-Electron Solid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.6 Reduced Density Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.7 Meaning of Single-Particle Eigenvalues in DFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.8 Exchange-Correlation Effects in the Wigner Crystal . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 4.9 Lindhard Dielectric Response Function for Free-Electron Gas . . . . . . . . . . . . . . . . . . . . . . . 117 4.10 Thomas-Fermi Screening Length at Zero-Temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 4.11 Pseudopotential Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.12 Ewald Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.13 Madelung Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 4.14 Ewald Summation of Madelung Energy for NaCl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 2
4.15 Car-Parrinello Lagrangian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5 Electronic Properties of Crystals 129 5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.2 Graphene Band Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 5.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 5.5 Electron and Hole Concentrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 5.6 Position of Fermi Level in Band Gap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 5.7 Semiconductor Doping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 5.8 More Realistic p-n Junction Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 5.9 Band Bending at Metal-Semiconductor Interfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6 Electronic Excitations 154 6.1 Plasma Frequency of Uniform Electron Gas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 6.2 Dielectric Function From Inter-Band Transitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 6.3 Dielectric Function From Intra-Band Transitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 6.4 Lorentz Model of Dielectric Function, Application to Si . . . . . . . . . . . . . . . . . . . . . . . . . . 157 6.5 Drude Model of Dielectric Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 6.6 Static Dielectric Constant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 6.7 Sum Rule for Imaginary Part of Dielectric Constant . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 6.8 Frenkel Exciton Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 6.9 Wannier Exciton Wavefunction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 6.10 Frenkel Excitons: Wavefunction and Bloch Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 6.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 7 Lattice Vibrations and Deformations 181 7.1 Phonons in Bulk Si Crystal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 7.2 Harmonic Oscillator Potential and Kinetic Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194 7.3 Phonon Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196 7.4 Thermal Expansion Coefficient, Grüneisen Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 198 7.5 Phonon Entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 7.6 UBER bulk modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202 7.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203 7.8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203 7.9 Elastic Constants of Isotropic Solid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203 7.10 Energy Density of Isotropic Solid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205 7.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206 7.12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206 7.13 Phonons of Graphene . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206 8 Phonon Interactions 219 8.1 Debye-Waller Factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219 8.2 Raman Scattering and Albrecht terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219 8.3 Cooper Pair Commutation Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221 8.4 BCS Number of Particles Variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221 8.5 BCS Gap Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222 9 Dynamics and Topological Constraints 226 9.1 Resistivity Tensor in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 226 9.2 Time-Reversal Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 226 9.3 Exponential Form of Time-Reversal Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 229 9.4 Berry Curvature of a Magnetic Monopole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230
3
9.5 Expectation Value For Electron Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230 9.6 Polarization of Finite System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232 9.7 Electron Velocity and Berry’s Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232 9.8 Berry Curvature of 2D Honeycomb Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233 9.9 Rice-Mele Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 234 9.10 T -Symmetry Breaking in Honeycomb Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 236 9.11 Dirac Points in Square Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 238 10 Magnetic Behavior of Solids 243 10.1 Hund’s Rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243 10.2 Magnetization and Band Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243 10.3 Stoner Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 245 10.4 Ising Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 248
4
Thanks to Laura Kulowski for undertaking a significant amount of compiling, equation number checking, and formatting of submitted solutions.
1
From Atoms to Solids
1.1
Packing Ratio
Solution by Robert Hoyt (2016) and Cedric Flamant (2018) An important consideration in the formation of crystals is the so-called packing ratio or filling fraction. For each of the elemental crystals (simple cubic, face-centered cubic, body-centered cubic), calculate the packing ratio, that is the percentage of volume occupied by the atoms modeled as touching hard spheres. Simple Cubic Note that the simple cubic structure, shown in Figure 1 (which can be tiled to cover all space) represents a unit volume taken up by the crystal. The cut spheres are the atoms. If the atom has radius r, then we see from the diagram that the side length of the unit block is 2r. Thus, the volume of the block, Vs , is 3
Vs = (2r) = 8r3 . The volume taken up by the atom can be determined by noting that there are 8 octants of the atom in the unit cell, thus giving the total volume of one atom. Thus, Va =
4 3 πr . 3
This leaves a packing ratio of the simple cubic structure ρsc of ρsc =
Va π = ≈ 0.5236. Vs 6
(1)
5
Figure 1: Simple Cubic
Face-Centered Cubic In the FCC structure we see that the radius of the unit cell is determined by the diagonal line of the face coinciding with four copies of the sphere’s radius. Thus, √ 3 Vs = 2 2r . There are also 8 octants and 6 hemispheres for a total of 4 spheres in this unit cell. Thus, there are a total of 4 spheres in the unit volume and Va =
16 3 πr . 3
This gives a packing ratio of π ρfcc = √ ≈ 0.7405. 3 2
(2)
Figure 2: Face-Centered Cubic
Body-Centered Cubic In the BCC structure we see that the diagonal of the cube is equal to four of the sphere’s radius. Hence, √ 3a2 = 4r r 16 2 4 a= r = √ r, 3 3 where a gives the cube side length, implying that 43 Vs = a3 = √ r 3 . 3 3 There are 8 octants and 1 full sphere in the unit cell for a total of 2 spheres. Thus, Va = This gives a packing ratio of
Figure 3: Body-Centered Cubic
8 3 πr . 3 √ π 3 ρbcc = ≈ 0.6802. 8
(3)
For the NaCl and CsCl structures, assuming that each type of ion is represented by a hard sphere and the nearest-neighbor spheres touch, calculate the packing ratio as a function of the ratio of the two ionic radii and find its extrema and the asymptotic values (for the ratio approaching infinity or zero, that is, with one of the two ions being negligible in size relative to the other one). Using values of the ionic radii for the different elements, estimate the filling fractions in the actual solids.
6
NaCl Structure The rock salt structure is technically two FCC lattices of Na and Cl shifted and superimposed on one another, as seen in Figure 4, but notice that simply an octant of this unit cell will suffice for analysis by symmetry. The first thing to notice about the rock salt structure is that there will be three regimes. Letting sodium be atom A and chloride be atom B and rA and rB their corresponding radii, we intuit that the contact points of the spheres will be different in the regions of rA rB , rA ∼ rB (in a manner yet to be made precise), and rA rB . Regime where rA rB In this regime, we see that the A atoms (blue) are much smaller than the B atoms (green). The size of this characteristic cell is determined by the Figure 4: Rock salt structure B atoms coming into contact with each other at the faces of the cube while the A atoms simply take up the free space in between the larger B atoms. We notice that the diagonal on√the face of the cube is coincident with 2 of √ √ 3 atom B’s radii, so this means 2a = 2rB , implying a = 2rB and Vs = a3 = 2 2rB . In total there are 4 octants of each atom in the cell, so Va = 4
14 3 2 14 3 3 3 πrA + 4 πrB = π rA + rB . 83 83 3
This gives a packing ratio of ρNaCl, =
Va π r3 + r3 π = √ A 3 B = √ 1 + x3 , Vs 3 2 rB 3 2
A where x ≡ rrB . Now we have to determine the limit of this regime. This occurs when atom A is large enough that it comes into contact with atom B along the faces as well—if it grows any bigger it will force atom B to shrink to make room. At this point, the side length of the cell will also be equal to rA + rB :
Figure 5: Regime where rA rB
√ rA = x = 2 − 1. rB √ Hence, this regime is described by x < 2 − 1. rA + rB =
√
2rB ⇒
Regime where rA ∼ rB In this regime, the size of the cell is determined by the contact of atom A and atom B along the edges of the cell. We see that the cell edge is equal to the sum of the two radii giving 3
Vs = (rA + rB ) . Hence we are left with a packing efficiency of 2 3 3 2 x3 + 1 Va 3 π rA + rB = ρNaCl,∼ = 3 = 3π 3. Vs (rA + rB ) (x + 1) Note that if we continue to increase the size of rA (blue) eventually we will reach the end of this regime—atom A will come into contact with its other selves. With similar reasoning to the last section, this happens when rA + rB =
√
2rA ⇒
7
rA 1 =x= √ , rB 2−1
which we could have predicted from the symmetry between atoms A and B. So, this regime exists when 1 . x ≤ √2−1
√
2−1 ≤
Regime where rA rB Now we consider the final regime where atom A has grown to dominate and hence determines the size of the cell. By symmetry, we recoginize that the packing efficiency must be Va π r3 + r3 π 1 ρNaCl, = = √ A 3 B = √ 1+ 3 . Vs x 3 2 rA 3 2 Now, putting the complete packing ratio as a function of ionic radii ratio together, we obtain √ π 3 x < 2 − 1, 3√2 1 + x √ 1 x3 +1 ρNaCl (x) = 2π , 2 − 1 ≤ x ≤ √2−1 3 3 (x+1) √ π √1 1+ 1 < x. 3 2
x3
Figure 6: Regime where rA ∼ rB
2−1
Plotting the above function, we get Figure 8, where the black dot shows the actual location of NaCl, using the ionic radii data found on the Wikipedia page of CsCl (it has a discussion of the CsCl compared to NaCl where the ionic radii at the relevant coordination number are mentioned) of xNaCl =
102 pm Na+ radius = . 181 pm Cl− radius
Explicitly, the estimated packing ratio of NaCl is ρNaCl = 0.646.
(4)
Figure 7: Regime where rA rB
It is worth mentioning that the packing ratio curve is understandable intuitively—consider Figure 9: In the first configuration, we see that the packing is not as efficient as it could be since the blue atoms (sodium) can be made bigger so that they completely fill the space between the chloride ions. Next, when the sodium atoms get too big to fill the holes, they start fighting for space with the chlorides and as such we see the formation of more gaps and a decrease in packing efficiency. Finally, when the sodium atoms are quite larger than the chlorides, we get a symmetric repeat of the first diagram. NaCl Packing Efficiency vs. Atomic Radii Ratio 0.80 0.75 0.70 0.65 0.60 0.55
1
2
3
4
Figure 8: NaCl Packing Efficiency Curve. The black dot indicates the actual packing efficiency of NaCl. The top green line shows the highest packing efficiency, and the lower orange line shows the asymptotic packing achieved when one atom is negligible in size relative to the other species.
8
Figure 9: NaCl structure with increasing ratio of atom A radius to atom B radius. π , and As for the points of interest of the packing ratio function: note the asymptotic values of ρNaCl (0) = 3√ 2 π ρNaCl (x → ∞) = 3√2 , which is the FCC packing ratio, as expected when one atom type is vanishingly small in two √ √ 1 interpenetrating FCC lattices; the maxima ρNaCl 2 − 1 = ρNaCl √2−1 = π 53 − 2 , which beats FCC packing.
The minimum is ρNaCl (1) = π6 , which corresponds to the Simple Cubic structure, which is what we get when both atom types are exactly the same size. CsCl Structure The CsCl structure is seen in Figure 10. Once again, we notice that there will be three regimes in the analysis of this structure. Letting cesium be atom A and chloride be atom B, we once again will have a region where atom A dominates, where atom B dominates, and a region where both play a role. Regime where rA rB Here we see that atom B determines the size of the cell by twice rB determining the edge length. This gives us 3 Vs = 8rB .
Figure 10: CsCl Structure
There are 8 octants of atom B and a full atom A in the cell, so Va =
4 3 3 rA + rB , 3
and we have a packing ratio of ρCsCl, =
π 3 x +1 . 6
The limit of this regime is when atom A grows large enough to touch the B atoms. This happens when the diagonal of the cube is also equal to twice the sum of the radii of the two atoms: √ √ rA 2 3rB = 2(rA + rB ) ⇒ = x = 3 − 1. rB
9
Figure 11: Regime where rA rB , rA ∼ rB , and rA rB , respectively. Regime where rA ∼ rB In this regime the diagonal of the cube is equal to twice the sum of the two ion’s radii, √
3a = 2(rA + rB ) 2 a = √ (rA + rB ), 3
8 3 ⇒ Vs = √ (rA + rB ) , 3 3
and hence √ π 3 x3 + 1 . ρCsCl,∼ = 2 (x + 1)3 1 This breaks down when x reaches √3−1 by symmetry with the previous situation—atom A will grow large enough that it touches the other atom A’s (not pictured, but it happens when A touches the bounding cell).
Regime where rA rB In this final regime, atom A is touching the bounding box, which means that it is pushing against the other atom A’s. By symmetry with the other extreme regime, π 1 ρCsCl, = + 1 . 6 x3 Putting it all together, we have the packing ratio function π 3 6 x + 1 √
ρCsCl (x) =
π 3 x3 +1
3
2 (x+1) π 1 + 1 6
x3
√ x ≤ 3 − 1, √ 1 3 − 1 < x < √3−1 , √1 ≤ x. 3−1
Plotting the above function, we get Figure 12, where the black dot shows the actual location of CsCl, xCsCl =
Cs+ radius 174 pm = . − 181 pm Cl radius
Explicitly, the estimated packing ratio of CsCl is ρCsCl = 0.6810.
(5)
Again we can understand the packing ratio function intuitively, as shown in Figure 13. In the first plot we can increase the size of the A atoms to fill up the free space. When the different atom types start competing with each 10
CsCl Packing Efficiency vs. Atomic Radii Ratio
0.70
0.65
0.60
0.55
1
2
3
4
Figure 12: CsCl Packing Efficiency Curve. Black dot shows the actual position of CsCl, and the green line at the top shows the maximum packing efficiency while the lower orange line shows the minimum asymptotic packing efficiency.
Figure 13: CsCl structures with increasing ratio of atom A radius to atom B radius. other for space, more gaps are formed, decreasing the efficiency. Then, when atom A becomes dominant, we have the symmetric behavior where atom B is relegated to filling up the gaps. As for the points of interest in the packing ratio function, note the asymptotic values of ρCsCl (0) = π6 , ρCsCl (x → ∞) = π NaCl and CsCl Packing and the maxima 0.80 6 , corresponding to Simple Cubic packing, √ √ 1 ρCsCl 3 − 1 = ρCsCl √3−1 = π 3 − 32 . The local mini0.75
√
mum is ρCsCl (1) = 83π , corresponding to BCC packing as the 0.70 atoms are the same size. Finally, before we move on to the next problem, consider the 0.65 NaCl and CsCl packing functions superimposed in Figure 14. On the CsCl page, Wikipedia says that when the ions are a 0.60 similar size, the CsCl structure is adopted, and when they are 0.55 of quite different sizes, the NaCl structure is used. The above 1 2 3 4 plot provides a possible justification for this statement. Note that the blue CsCl curve has a higher packing efficiency for ionic radii ratios near unity, while the red NaCl curve does better away Figure 14: NaCl and CsCl packing efficiency from one. Note that the locations of the real salts are on the curves superimposed. appropriate curves to maximize their own packing efficiencies. Neat.
11
1.2
Madelung Energy of 2D Bipartite Lattices
Solution by Cedric Flamant (2018) The three ionic lattices, rock salt, cesium chloride and zinc blende, are called bipartite lattices because they include two equivalent sites per unit cell which can be occupied by the different ions so that each ion type is completely surrounded by the other. Describe the corresponding bipartite lattices in 2 dimensions. Are they all different from each other? Using the interactive structures shown on http://www.chemtube3d.com/solidstate/_table.htm, we can find the 2D projected bipartite lattices. These are shown in Figure 15.
Figure 15: 2D bipartite lattice projections. From left to right, NaCl, CsCl, Zinc blende 1 and 2. All three share a common 2D bipartite lattice, but zinc blende has an additional honeycomb lattice. Try to obtain the Madelung energy for one of them, and show how the calculation is sensitive to the way in which the infinite sum is truncated. The Madelung energy is the energy associated with picking a specific ion in the lattice and computing its interaction energy with all the other ions in the lattice. Formally for the square AB lattice, it is given by ME =
0 ∞ X q2 (−1)i+j 4π 0 r0 i,j=−∞ (i2 + j 2 )1/2
where the prime on the sum indicates i = j = 0 is to be left out. The r0 designates the nearest-neighbor distance. The Madelung constant is purely geometrical and is the result of the infinite sum. We can see that the summand’s numerator makes sense since it gets the attraction/repulsion right depending on which ion pair is being considered. Now the question is how the calculation is sensitive to summing order. Note that this is what is meant by “sensitive to the way in which the infinite sum is truncated” since in general we cannot add every term in the series, we will have to stop at some point. Which terms we include matter, and we will see that some sums head towards convergence nicely, while others do not. For the sake of this question, only showing how one of the converging sums depends on truncation is insufficient, since this does not exhibit the sensitivity of this mathematical expression. There are two primary candidates that first come to mind: summing by squares and summing by circles. The two approaches to summing the energy contribution are shown in Figure 16. Squares: The idea here is to sum outwards, adding up all new charges that fall within squares of increasing size. Formally, this would be represented with the mathematical expression lim
r→∞
0 r i+j X (−1) 2 2 1/2 i,j=−r (i + j )
12
.
Figure 16: Expanding by circles and expanding by squares on a square lattice. Circles: Here the idea is to sum outwards, adding up all charges that fall in successively larger circles. Formally, this would correspond to the mathematical expression (which can be directly translated into code), r (−1)i+j , i2 + j 2 1/2 ≤ r X (i2 +j 2 )1/2 lim . 1/2 r→∞ 0, i = j = 0 or i2 + j 2 >r i,j=−r
These two approaches give wildly different convergence behavior, as seen in the following Mathematica notebook, Figure 17. Note the very different convergence behavior of the two methods. The discrepancy can be understood in terms of the balance of charge added in each increase of the size of the expanding square/circle. For the square, a neutral total charge is added each time more ions are considered, so as you go farther out the terms are smaller corrections. Only the Coulomb 1/r factor is relevant, and hence the additional terms decreasing in size. But, for the circle method, the charge imbalance per increase in radius tends to grow, in fact, it is essentially growing as fast as the 1/r fall-off of the potential! Note that in Figure 16, in successive circles the number of charges included goes up roughly as the circumference of the circle, 2πr, and they will always be predominantly positive, or negative. So, additional terms are (to leading order) proportional to r/r = const.! That would explain the very slow and haphazard approach to the Madelung constant. In fact, it is not evident that it is a convergent series at all—in the 3D equivalent, NaCl, expanding spheres has been proven to not converge. Other summing orders exist, and some will converge nicely while others will not.
13
Figure 17: Plot showing the convergence difference between expanding squares and circles. In blue we have the smooth convergence of expanding squares, and the orange dots are expanding circle calculations as a function of radius.
1.3
Simultaneous Energy and Momentum Eigenstates
Solution by Daniel Larson (2019) Consider the single-particle hamiltonian: Ĥsp =
p2 + V(r) 2m
where p is the momentum operator p = −i~∇r m the mass of the particle, and V(r) the potential energy. Show that the commutator of the hamiltonian with
14
the potential energy term, defined as h
i Ĥsp , V(r) = Ĥsp V(r) − V(r)Ĥsp
is not zero, except when the potential is a constant V(r) = V0 . Based on this result, provide an argument to the effect that the energy eigenfunctions can be simultaneous eigenfunctions of the momentum operator only for free particles. We begin by carefully expanding and computing the various terms in the commutator. 2 2 h i p2 p p sp Ĥ , V(r) = + V(r), V(r) = , V(r) + [V(r), V(r)] = , V(r) 2m 2m 2m because the potential commutes with itself. The commutator of a differential operator like p is defined by its action on an arbitrary test function f (r).
2 p2 −~ 2 −~2 2 −~2 2 , V(r) f (r) = ∇r , V(r) f (r) = ∇r (V(r)f (r)) − V(r) ∇ (f (r)) 2m 2m 2m 2m r
(6)
The first term on the right hand side of (6) can be expanded using vector derivative identies: ∇2r (V(r)f (r)) = ∇2r V f + 2 (∇r V) · (∇r f ) + V∇2r f The last term above cancels the final term in (6). Suppressing the arbitrary test function f , we arrive at the final result for the commutator: h i −~2 ∇2r V(r) + 2 (∇r V(r)) · ∇r Ĥsp , V(r) = 2m This commutator will only vanish for all r if V is a constant V0 . Because a vanishing commutator is a requirement for two operators to have simultaneous eigenfunctions, energy eigenfunctions of the single-particle hamiltonian will only be simultaneous eigenfunctions of the momentum operator if h i 0 = Ĥsp , p = [V(r), p] = i~ (∇r V(r)) , where we have used the same procedure as above to evaluate this second commutator. It will only vanish if V(r) = V0 is a constant. In that case, the single-particle hamiltonian represents a free particle with the energy of a stationary particle (p = 0) given by V0 .
1.4
Bond Formation in 1D
Solution by Robert Hoyt (2016) In the example of bond formation in a 1D molecule, using the definitions of the wavefunctions for the isolated atoms, Eq. (1.14), and the bonding and anti-bonding states, Eq. (1.15), show that: (a) The difference in the probability of an electron being in the bonding or anti-bonding states, rather than in the average of the two isolated-atom wavefunctions, is given by the expression in Eq. (1.16).
15
b 1 ψ1 (x) = √ e−|x− 2 |/a a 1 −|x+ b |/a 2 ψ2 (x) = √ e a b b 1 ψ ± (x) = √ [e−|x− 2 |/a ± e−|x+ 2 |/a ] N± 1 (ψ1 (x) ± ψ2 (x)) =q 2(1 ± e−b/a (1 + ab ))
(7a) (7b)
(7c)
The probability of the electron in the symmetric or the antisymmetric state is: (ψ ± (x))∗ ψ ± (x) =
1 2(1 ± e−b/a (1 + ab ))
|ψ1 (x)|2 + |ψ2 (x)|2 ± (ψ1 (x)ψ2 (x) + c.c.)
(8)
The probability to be in the average of the states ψ1 (x) and ψ2 (x) is : 1 (|ψ1 (x)|2 + |ψ2 (x)|2 ) 2
(9)
and the difference is: b 2 2 −b/a )(|ψ (x)| + |ψ (x)| ) ±2ψ (x)ψ (x) ∓ e (1 + 1 2 1 2 a 2 1 ± e−b/a 1 + ab a b = ± 2ψ1 (x)ψ2 (x) − e−b/a (1 + ) |ψ1 (x)|2 + |ψ2 (x)|2 λ a 1
δn± =
(10)
where we have used the fact that the wavefunctions are real. (b) The change in the potential energy for the combined system versus the two isolated ions, given by Eq. (1.17), is negative for the ψ (+) state and positive for the ψ (−) state (take b = 2.5 a and do a numerical integration). (c) The change in the kinetic energy for the combined system versus the two isolated ions is negative for the ψ (+) state and positive for the ψ (−) state (take b = 2.5 a and do a numerical integration).
1 ∆V = ψ V1 (x) + V2 (x) ψ − hψ1 |V1 (x)|ψ1 i + hψ2 |V2 (x)|ψ2 i (11) 2 1 (12) = ± hψ1 | ± hψ2 | (V1 (x) + V2 (x)) |ψ1 i ± |ψ2 i − hψ1 |V1 (x)|ψ1 i N 1 = ± hψ1 |V1 (x) + V2 (x)|ψ1 i ± hψ1 |V1 (x) + V2 (x)|ψ2 i + hψ2 |V1 (x) + V2 (x)|ψ2 i ± hψ2 |V1 (x) + V2 (x)|ψ1 i N (13) ±
±
− hψ1 |V1 (x)|ψ1 i 1 ∆V = ± 2 − N ± hψ1 |V1 (x)|ψ1 i + 2 hψ1 |V2 (x)|ψ1 i ± 4 hψ1 |V1 (x)|ψ1 i . N
(14) (15)
In Eq. (12) through (15), the fact that the potentials are real and the symmetry about the origin x = 0 was used to greatly simplify the problem. See the following Mathematica code for the numerical integration.
16
Here are the definitions I used, including gaussian pseudopotentials. a = 1; b = 25 /10; NN[s_] := 2 1 +s *Exp-b a 1 +b a; 1 Exp-Absx -b 2a; p1[x_] := Sqrt[a] 1 Exp-Absx +b 2a; p2[x_] := Sqrt[a] -1 2 Exp-2 x -b 2 ; V1 = Sqrt[π/2] -1 2 Exp-2 x +b 2 ; V2 = Sqrt[π/2]
Normalization constants for the symmetric (NN[1]) and antisymmetric (NN[-1]) states: NN[1] // N NN[-1] // N 2.57459 1.42541
Change in the potential energy: First term, <1|V1|1>, the “on-site” integral onsite= NIntegrate[p1[x]*V1 *p1[x], {x, -10, 10}] -0.523157
Second term, <1|V2|1>, the “off-site” integral. It’s much smaller. offsite= NIntegrate[p1[x]*V2 *p1[x], {x, -10, 10}] -0.0111089
Third term, <1|V1|2>, the “diagonal” integral. It’s in the middle. diagonal= NIntegrate[p1[x]*V1 *p2[x], {x, -10, 10}] -0.0625141
Putting it all together. Since the equation is non-dimensionalized, there is an implicit factor of e2 (4 π ϵ0 a) = 2 Ry scaling these values (MKS units). ΔV+ = ΔV- =
1 NN[1] 1
(2 -NN[1])*onsite+2 *offsite+4 *diagonal
NN[-1]
(2 -NN[2])*onsite+2 *offsite-4 *diagonal
0.0110032 0.581621
It’s important to note that if the atoms were farther apart, the normalization constant approaches 2 and the expression simplifies to (offsite ± 2*diagonal). Since the diagonal term is then always larger than the offsite term, the potential change becomes strictly negative for the symmetric state.
Printed by Wolfram Mathematica Student Edition
17
2
1D_hydrogen_numerical.nb
potential change
strictly negative
symmetric
When the atoms are somewhat close together like they are here, however, the state must remain normalized so some density is removed from the region where the potential is most negative. This yields the slightly positive potential energy change. This might appear to disprove the molecular bonding hypothesis, but we’ve neglected one very important fact: the wavefunctions did not distort at all. In more sophisticated calculations we would variatonally optimize the shape of the wavefunctions (i.e. allow them to distort) in addition to forming symmetric and antisymmetric combinations, yielding lowerenergy solutions.
Change in kinetic energy There are two ways to approach this problem. The first is the following: Mathematica doesn’t directly differentiate Abs[x ± b/2a] symbolically, so I use a combination of Assuming and Simplify to allow Mathematica to eliminate the Abs[...] factors in the wavefunctions before calculating the gradient. At the end, I define a piecewise function to stitch the three solution regions together. The steps go something like this under each Assuming command: 1. simplify the wavefunctions, eliminating the Abs factors 2. take the second derivative of the Abs-free wavefunctions 3. multiply by another factor of the wavefunctions (the bra) 4. simplify the overall expression to get rid of the Abs[...] factors introduced in step 3 Everything is left as a function of s, where s=+1 is the symmetric state and s=-1 is the antisymmetric state. KElow[x_, s_] := Assuming u < -b 2 a,
-1 NN[s]
(p1[u]+s p2[u])*D(p1[u]+s p2[u]) // Simplify , {u, 2}+
1
p1[u] Dp1[u] // Simplify , {u, 2}+p2[u] Dp2[u] // Simplify , {u, 2} // 2 Simplify /. u → x KEmid [x_, s_] := Assuming u > -b 2 a && u < b 2 a, -1 NN[s]
(p1[u]+s p2[u])*D(p1[u]+s p2[u]) // Simplify , {u, 2}+ 1
p1[u] Dp1[u] // Simplify , {u, 2}+p2[u] Dp2[u] // Simplify , {u, 2} // 2 Simplify /. u → x KEhigh[x_, s_] := Assuming u > b 2 a, -1 NN[s]
(p1[u]+s p2[u])*D(p1[u]+s p2[u]) // Simplify , {u, 2}+ 1
p1[u] Dp1[u] // Simplify , {u, 2}+p2[u] Dp2[u] // Simplify , {u, 2} // 2 Simplify /. u → x KEtot[x_, s_] := PiecewiseKElow[x, s], x < -b 2 a, KEmid [x, s], x > -b 2 a && x < b 2 a, KEhigh[x, s], x > b 2 a
Printed by Wolfram Mathematica Student Edition
18
1D_hydrogen_numerical.nb
3
Plot{KEtot[x, -1], KEtot[x, 1]}, {x, -5, 5}, PlotRange → Full, ImageSize → Small 0.05
-4
2
-2
4
-0.05
The expected change in kinetic energy turns out to be zero for both wavefunctions! The values are just numerical error, shown by the factor 10-90 introduced by integrating more accurately. NIntegrateKEtot[x, +1], {x, -10000, +10 000}, WorkingPrecision→ 80, AccuracyGoal → 80, PrecisionGoal→ 80 NIntegrate[KEtot[x, -1], {x, -10, +10}] -1.6303426025503643046633099954781661831787430873110051479875163112017898047973565 ×10-99 2.20312 ×10-9
The second approach is to fight smarter rather than harder. Taking the second derivative of the ψi ’s directly reveals an interesting fact about the single-atom wavefunctions... In[75]:=
ψ1 [x_] := 1
PiecewiseExp-x -b 2a, x > b 2 a, Expx -b 2a, x ≤ b 2 a; Sqrt[a] ψ1 [x] ψ1 ''[x] ψ1 ''[x]-ψ1 [x] // Simplify 5
x > 54
5 - 4 +x
x ≤ 54 True
ⅇ 4 -x Out[76]=
Out[77]=
Out[78]=
ⅇ 0
ⅇ- 4 +x
5
x < 54
5 -x 4
ⅇ Indeterminate
x > 54 True
Indeterminate 0
4x ⩵5 True
Aside from the discontinuity at x=b/2a, the isolated atom wavefunctions are eigenstates of the kinetic energy operator up to a factor of 1a2 introduced by the second derivative. The symmetric/antisymmetric wavefunctions are linear combinations of these wavefunctions, and are therefore also eigenstates. 1
1
1
Thus the kinetic energy difference reduces to a2 times (norm of ψ±) - 2 (norm of ψ1 ) - 2 (norm of ψ2 ) = 1 1 1 1 - 2 - 2 = 0, i.e. zero change in the kinetic energy. a2
Again, the explanation for this apparent problem is that the isolated-atom wavefunctions remained unchanged. When the wavefunctions are variationally optimized, the distortion from the isolated wavefunctions produces a significant change in the kinetic energy.
Printed by Wolfram Mathematica Student Edition
19
1.5 1.6
Bonding in Water
Solution by Zoe Zhu (2019) Produce an energy-level diagram for the orbitals involved in the formation of the covalent bonds in the water molecule. Provide an argument of how different combinations of orbitals than the ones discussed in the text would not produce as favorable a covalent bond between H and O. Describe how the different orbitals combine to form hydrogen bonds in the solid structure of ice.
Figure 18: Energy levels of H2 O. The s orbital comes from the hydrogen atom and the p orbital comes from the oxygen atom, so the hybrid bond must be an spn -type hybrid. However, an sp hybrid, meaning one s orbital mixes with one of the three p orbitals, is linear. Recall that there are two covalent bonds, and therefore we cannot have a linear configuration between the bonds. Therefore, only sp2 and sp3 hybrids are possible. In water, hydrogen bonds are constantly breaking and reforming. In ice, however, the crystalline structure is maintained by hydrogen bonding because there is no longer enough energy to break the hydrogen bond. In ice, H2 O molecules are orientational and the hydrogen bonds are longer due to its weak nature, causing ice to be denser than its liquid form. Ice can have different crystal structures. In the cubic ice, each oxygen atom is connected to four hydrogen atoms by two covalent bonds and two hydrogen bonds, forming an sp3 hybrid. In the hexagonal ice, each oxygen atom is connected to three hydrogen atoms by two covalent bonds and one hydrogen bonds, forming an sp2 hybrid.
1.7
Excitations of the Free-Electron Model
Solution by Robert Hoyt (2016)
20
Consider a simple excitation of the ground state of the free-electron system, consisting of taking an electron from a state with momentum k1 and putting it in a state with momentum k2 ; since the ground state of the system consists of filled single-particle states with momentum up to the Fermi momentum kF , we must have |k1 | ≤ kF and |k2 | > kF . Removing the electron from state k1 leaves a “hole” in the Fermi sphere, so this excitation is described as an “electron-hole pair”. Discuss the relationship between the total excitation energy and the total momentum of the electron-hole pair; show a graph of this relationship in terms of reduced variables, that is, the excitation energy and momentum in units of the Fermi energy F and the Fermi momentum kF . (At this point we are not concerned with the nature of the physical process that can create such an excitation and with how momentum is conserved in this process.) First, note that in a full Fermi sphere there is no net momentum due to spherical symmetry, ksphere = 0. However, imagine a hole where k1 would be. Then, since kholesphere + k1 = ksphere = 0, we deduce that kholesphere = −k1 . Adding on the new momentum of the excited electron gives a net momentum k of the electron-hole pair of k = k2 − k1 . As for the excitation energy, it will simply be the energy difference between the new state and the old (hole) state: E=
~2 2 2 |k2 | − |k1 | . 2m
For cleaner notation, and because the problem ultimately wants the relationship in reduced units anyways, we can re-express E, k, k1 and k2 as E/ F → E Noting that F =
and k/kF → k.
2 ~2 k F 2m , our earlier equations become
k = k2 − k1
(16)
2
(17)
k = k22 + k12 − 2k1 · k2 E = k22 − k12 .
(18)
Substituting the identity k2 = k + k1 into the expression for E we get E = k 2 + 2k · k1 . For clarity in discussion, we define cos θ as the vector angle between k and k1 to obtain E = k 2 + 2kk1 cos θ. The range of energy values depends on k1 and cos θ. Since k1 ≥ 0 is the vector magnitude of k1 , extrema occur for k1 = 1 where the hole state is right at the Fermi level. This simplification yields E = k 2 + 2k cos θ. The maximum range is controlled by cos θ, whose extreme values are ±1. For sufficiently small values of k we would have k 2 < 2k, which for negative values of cos θ would lead to a negative energy E. However, thinking about the Fermi sphere, we note that small k would necessarily correspond to an electron right under the surface being excited to a state right outside the surface of the sphere. Thinking back to the definition of θ, the angle between k and k1 , a negative value of cos θ implies that the electron moves deeper into the Fermi sphere. This is not allowed since the interior of the Fermi sphere is already filled. However, we see that we could have an arbitrarily small excitation to a state outside the sphere. As k gets larger though, we see that E can only remain infinitesimally close to zero until k = 2. Beyond this, E must become nonzero. This limit can be
21
visualized as an excitation that just barely traverses the Fermi sphere to the antipode, since the Fermi sphere has a diameter of 2 in reduced units. The rest of the boundaries are straightforward to determine using the full range of cos θ, and are given by
Emin ≈ 0
k≤2
(19)
Emin = k(k − 2)
k>2
(20)
Emax = k(k + 2)
∀k.
(21)
The allowed region is shaded in the figure below. E 35 30 25 20 15 10 5 1
1.8
2
3
4
5
k
Attractive Component of Lennard-Jones
Solution by Cedric Flamant (2018) In order to derive the attractive part of the Lennard-Jones potential, we consider two atoms with Z electrons each and filled electronic shells. In the ground state, the atoms will have spherical electronic charge distributions and, when sufficiently far from each other, they will not interact. When they are brought closer together, the two electronic charge distributions will be polarized because each will feel the effect of the ions and electrons of the other. We are assuming that the two atoms are still far enough from each other so that their electronic charge distributions do not overlap, and therefore we can neglect exchange of electrons between them. Thus, it is the polarization that gives rise to an attractive potential; for this reason this interaction is sometimes also referred to as the “fluctuating dipole interaction”. To model the polarization effect, we consider the interaction potential between the two neutral atoms: Vint =
X Ze2 Z 2 e2 − (1) |R1 − R2 | r −R i
−
2
i
X
Ze2
j
rj − R1
+
(2)
(1)
e2
X (1)
ri
ij
(2)
− rj
(2)
where R1 , R2 are the positions of the two nuclei and ri , rj are the sets of electronic coordinates associated with each nucleus. In the above equations, the first term is the repulsion between the two nuclei, the second term is the attraction of the electrons of the first atom to the nucleus of the second, the third term is the attraciton of the electrons of the second atom to the nucleus of the first, and the last term is the repulsion between the two sets of electrons in the two different atoms. From second order perturbation theory, the energy change due to this interaction is given by: D ∆E =
D
(1)
(2)
Ψ0 Ψ0
(1)
(2)
Vint Ψ0 Ψ0
E
+
X nm
(1)
(2)
(1)
(2)
Ψ0 Ψ0
(1)
(2)
Vint Ψn Ψm
E2
E0 − Enm (1)
(2)
where Ψ0 , Ψ0 are the ground-state many-body wavefunctions of the two atoms, Ψn , Ψm are their excited states, and E0 , Enm are the corresponding energies of the two-atom system in their unperturbed states.
22
We define the electronic charge density associated with the ground state of each atom through: Z 2 (I) (I) Ψ0 (r, r2 , r3 , . . . , rZ ) dr2 dr3 · · · drZ n0 (r) = Z =
Z Z X
2
(I)
δ(r − ri ) Ψ0 (r1 , r2 , . . . , rZ ) dr1 dr2 · · · drZ
i=1
with I = 1, 2 (the expression for the density n(r) in terms of the many-body wavefunction |Ψi is discussed in detail in Appendix A). Show that the first order term in ∆E corresponds to the electrostatic interaction energy between these two charge densities, show that this term vanishes (the two charge densities in the unperturbed ground state are spherically symmetric). We are asked to compute the first order energy perturbation, D E (1) (2) (1) (2) Ψ0 Ψ0 Vint Ψ0 Ψ0 . Let’s do it term by term in Vint . For the first term, D
(1)
(2)
Ψ0 Ψ0
E Z 2 e2 Z 2 e2 (1) (2) Ψ0 Ψ0 = , |R1 − R2 | |R1 − R2 |
where the operator is actually just a scalar since it only depends on the fixed positions of the nuclei. For the second term, D X (1) (2) Ψ0 Ψ0 − i
Ze2 (1) |ri − R2 |
(1)
(2)
Ψ0 Ψ0
E
= −Ze2
Z Z X
2
(1)
Ψ0 (r1 , . . . , rZ )
i=1
= −Ze2
Z Z X
1 dr1 · · · drZ |ri − R2 | 2
(1)
δ(r − ri ) Ψ0 (r1 , . . . , rZ )
i=1 2
Z
= −Ze
= −Ze2
1 dr dr1 · · · drZ |r − R2 |
Z Z X 2 1 (1) dr δ(r − ri ) Ψ0 (r1 , . . . , rZ ) dr1 · · · drZ |r − R2 | i=1 (1)
Z dr
n0 (r) . |r − R2 |
Similarly, for the third term, D X (1) (2) Ψ0 Ψ0 − j
Ze2 (2)
|rj − R1 |
(1)
(2)
Ψ0 Ψ0
23
E
= −Ze2
(2)
Z dr
n0 (r) . |r − R1 |
Finally for the last electron-electron interaction term, D
(1)
(2)
Ψ0 Ψ0
e2
X
(1)
(2)
(1)
(2)
Ψ0 Ψ0
E
|ri − rj | Z 2 2 X (2) (2) (2) (1) (1) (1) Ψ0 r1 , . . . , rZ Ψ0 r1 , . . . , rZ = e2 ij
ij
= e2
XZ
(1)
Ψ0
(1)
(1)
r1 , . . . , rZ
i
= e2
Z Z 2X
(1)
Ψ0
Z dr1
Z Z X i=1
= e2
(2)
dr2 δ r2 − ri
(1)
(1)
(1)
r1 , . . . , rZ
2Z
Z dr1
dr2
(2)
(2)
(1) r1 − r2
(1)
(1)
(2)
(2)
dr1 · · · drZ dr1 · · · drZ
(2)
dr2
n0 (r2 ) (1) r1 − r2
(1)
(1)
dr1 · · · drZ
Z (2) 2 n (r2 ) (1) (1) (1) (1) (1) (1) Ψ0 r1 , . . . , rZ dr1 · · · drZ dr2 0 δ r1 − ri (1) r1 − r2 (1)
Z
(1)
dr1 · · · drZ dr1 · · · drZ
2 (2) (2) Ψ(2) r1 , . . . , rZ 0
j=1
XZ i
= e2
1 (1) (2) ri − rj
(2)
n0 (r1 )n0 (r2 ) . |r1 − r2 |
So, in total, Z 2 e2 ∆E = − Ze2 |R1 − R2 | (1)
Z
(1)
n (r) dr 0 − Ze2 |r − R2 |
Z
(2)
n (r) dr 0 + e2 |r − R1 |
Z
(1)
Z dr1
dr2
(2)
n0 (r1 )n0 (r2 ) . |r1 − r2 |
(2)
The spherical symmetry of n0 (r) and n0 (r) allows us to apply Gauss’s law, which tells us that the terms are equivalent to point charges, with charge equal to the integral of the total electron charge densities, centered at the position of each nucleus. This leaves us with Z 2 e2 Z 2 e2 Z 2 e2 ∆E = − − + e2 |R1 − R2 | |R1 − R2 | |R1 − R2 | Z (1) Z 2 e2 Zn0 (r1 ) =− + e2 dr1 |R1 − R2 | |r1 − R2 | 2 2 2 2 Z e Z e + =− |R1 − R2 | |R1 − R2 |
Z
(1)
Z dr1
dr2
(2)
n0 (r1 )n0 (r2 ) |r1 − r2 |
= 0, where again in the second and third equality we have made use of the equivalency granted by Gauss’s law given that the charge distributions are spherically symmetric about their respective nuclei, and that the charge distributions of the two atoms don’t overlap. We have shown that the first order term of the interaction vanishes, so we have to proceed to second order. The wavefunctions involved in the second order term in ∆E will be negligible, unless the electronic coordinates associated with each atom are within the range of non-vanishing charge density. This implies that the distances (1) (2) ri − R1 and rj − R2 should be small compared to the distance between the atoms |R2 − R1 |, which defines the distance at which interactions between the two charge densities becomes negligible. Show that (1) (2) expanding the interaction potential in the small quantities ri − R1 /|R2 − R1 | and rj − R2 /|R2 − R1 |,
24
gives, to lowest order: (2) r − R · (R2 − R1 ) − R · (R − R ) 2 1 2 1 j e · Vint ≈ − 3 2 2 |R2 − R1 | ij (R2 − R1 ) (R2 − R1 ) (1) (2) X r i − R1 · r j − R2 e2 + . 2 |R2 − R1 | ij (R2 − R1 )
2
X
(1)
ri
Credit to Robert Hoyt (2016) for this section. So, we have to go to the second term of ∆E, which is second order perturbetion. To get the leading term, we first r
(1)
r
−R
(2)
−R2
R2 −R1 expand Vint in xi = |Ri 2 −R11| , yj = |Rj 2 −R1 | and R = |R1 − R2 |, R̂ = |R . 1 −R2 |
Vint =
Z Z Z 1 X Ze2 1 X Ze2 1 X Ze2 e2 − + − R R i=1 |xi − R̂| R i=1 |yi + R̂| R i=1,j=1 |xi − yj − R̂|
(22)
We can use the following vector Taylor expansion, which can be straightforwardly derived from the standard multidimensional Taylor expansion in Cartesian coordinates: 1 1 ~ 1 + 1 (r · ∇) ~ 2 1 + ... = ± (r · ∇) |R ± r| R R 2 R 2 1 r · R 3(r · R) r2 = ∓ + − + .... 3 5 R R 2R 2R3
(23)
Using this expansion, to second order in the small terms xi and yi , Vint =
Z 1 i Z 2 e2 1 X 2h 3 − Ze 1 + xi · R̂ + (R̂ · xi )2 − x2i R R i=1 2 2 Z
−
3 1 i 1 X 2h Ze 1 − yi · R̂ + (R̂ · yi )2 − yi2 R i=1 2 2
+
i 1 X 2h 3 1 Ze 1 + R̂ · (xi − yj ) + (R̂ · (xi − yj ))2 − (xi − yj )2 + O x3 R i,j 2 2
Z
(24)
0th and 1st order terms in the small distances are all cancelled and we are left with the second order terms:
Vint ≈ −
Z
Z
Z
Z
1 X X 2 3 3 3 1 1 1 e (R̂ · xi )2 + (R̂ · yi )2 − (R̂ · xi − R̂ · yj )2 − x2i − yi2 + (xi − yj )2 R i=1 j=1 2 2 2 2 2 2
i e2 X X h =− 3(R̂ · xi )(R̂ · yj ) − xi · yj R i=1 j=1 (1) (2) (1) (2) Z X Z h X (ri − R1 ) · (R2 − R1 )(rj − R2 ) · (R2 − R1 ) (ri − R1 ) · (rj − R2 ) i e2 =− 3 − |R2 − R1 | i=1 j=1 |R2 − R1 |4 |R2 − R1 |2
(25)
Using the above approximation for Vint , show that the leading order term in the energy difference ∆E behaves like −6 |R2 − R1 | and is negative. This establishes the origin of the attractive term in the Lennard-Jones potential. 1 So, now note that the interaction potential is proportional to |R −R 3 . Recalling the second order perturbation 2 1| to energy,
∆E ≈
D E2 1 (1) (2) (2) Ψ0 Ψ0 Vint Ψ(1) Ψ , n m E0 − Enm nm
X
25
E (1) (2) we note that in general Ψn Ψm will not exhibit the spherical symmetry around each nucleus that caused terms E D (1) (2) (1) (2) to cancel when evaluating Ψ0 Ψ0 Vint Ψ0 Ψ0 , so the sum in the second order perturbation is expected to be 1 nonzero. We then know that the full correction will be proportional to |R −R due to the squaring of the matrix |6 2
1
1 element in the expression. Finally, since the excited states in general have energy Enm > E0 , E0 −E < 0 and hence nm the above energy term is also negative. This establishes the attractive ∝ R16 term of the Lennard-Jones potential.
1.9
Approximate Morse Potential Eigenvalues
Solution by Cedric Flamant (2019) We wish to determine the eigenvalues of the Morse potential, Eq. (1.19). One method is to consider an expansion in powers of (r − r0 ) near the minimum and relate it to the harmonic oscillator potential with higher-order terms. Specifically, the potential V(r) =
1 mω 2 (r − r0 )2 − α(r − r0 )3 + β(r − r0 )4 2
(26)
1 1 ~ω 1 − γ n + n = n + 2 2
(27)
has eigenvalues
where 3 γ= 2~ω
~ mω
2
5 α2 −β . 2 mω 2
(28)
First, check to what extent the expansion (26) with up to fourth-order terms in (r − r0 ) is a good representation of the Morse potential; what are the values of α and β in terms of the parameters of the Morse potential? Use this approach to show that the eigenvalues of the Morse potential are given by Eq. (1.22). Starting the with Morse potential and expanding about the minimum at r0 , h i VM (r) = e−2(r−r0 )/b − 2e(r−r0 )/b (r − r0 )2 (r − r0 )3 7(r − r0 )4 − + ≈ −1 + b2 b3 12b4 (r − r0 )3 7 (r − r0 )4 (r − r0 )2 − + . = − + 2 3 b b 12b4
(29) (30) (31)
As long as |r − r0 | < b it is a good approximation, as can be verified by plotting. By comparison with Eq. (26), we 7 1 2 see that α = b 3 and β = 12 b4 , and b2 = 2 mω . Plugging these values in to the expression for γ, we find γ=
~ω , 4
which when inserted to the expression for the eigenvalues, ~ω 1 1 ~ω 1 − n+ , n = n + 2 4 2 which are the approximate eigenvalues of the Morse potential.
26
(32)
(33)
1.10
Kronig-Penney Model
Solution by Robert Hoyt (2016) An important simple model that demonstrates some of the properties of electron states in infinite-periodic solids is the so-called Kronig-Penney model. In this model, a particle of mass m experiences a one-dimensional periodic potential with period a: V(x) = 0,
0 < x < (a − l)
(34)
= V0 ,
(a − l) < x < a
(35)
V(x + a) = V(x)
(36)
where we will take V0 > 0. The wavefunction ψ(x) obeys the Schrödinger equation ~2 ∂ 2 + V(x) ψ(x) = ψ(x). − 2m ∂x2
(37)
(a) Choose the following expression for the particle wavefunction: ψ(x) = eikx u(x)
(38)
and show that the function u(x) must obey the equation ∂ 2 u(x) 2m 2mV(x) ∂ 2 u(x) 2 + 2ik − k − 2 + u(x) = 0. ∂x2 ∂x2 ~ ~2
(39)
Assuming that u(x) is finite for x → ±∞, the variable k must be real so that the wavefunction ψ(x) is finite for all x. (b) We first examine the case > V0 > 0. Consider two solutions u1 (x), u2 (x) for the ranges 0 < x < (a − l) and (a − l) < x < a, respectively, which obey the equations ∂ 2 u1 (x) ∂u1 (x) 2 + 2ik − k − κ2 u1 (x) = 0, 2 ∂x ∂x ∂u2 (x) 2 ∂ 2 u2 (x) + 2ik − k − λ2 u2 (x) = 0, 2 ∂x ∂x
0 < x < (a − l)
(40)
(a − l) < x < a
(41) (42)
where we have defined the quantities r κ=
2m , ~2
r λ=
2m( − V0 ) ~2
(43)
which are both real for > 0 and ( − V0 ) > 0. Show that the solutions to these equations can be written as: u1 (x) = c1 ei(κ−k)x + d1 e−i(κ+k)x ,
0 < x < (a − l)
(44)
−i(λ+κ)x
(a − l) < x < a.
(45)
u2 (x) = c2 e
i(λ−κ)x
+ d2 e
,
By matching the values of these solutions and of their first derivatives at x = 0 and x = a − l, find a system of four equations for the four unknows, c1 , d1 , c2 , d2 . Show that requiring this system to have a non-trivial solution leads to the following condition: −
κ2 + λ2 sin (κ(a − l)) sin (λl) + cos (κ(a − l)) cos (λl) = cos(ka). 2κλ
27
(46)
Next, show that with the definition 2 κ + λ2 tan(θ) = − tan (λl) 2κλ
(47)
the above condition can be written as: "
κ2 − λ2 1+ 4κ2 λ2
#1/2
2 2
sin (λl)
cos (κ(a − l) − θ) = cos (ka).
(48)
Show that this last equation admits real solutions for k only in certain intervals of ; determine these intervals of and the corresponding values of k. Plot the values of as a function of k and interpret the physical meaning of these solutions. How does the solution depend on the ratio /V0 ? How does it depend on the ratio l/a? (c) Repeat the above problem for the case V0 > > 0. Discuss the differences between the two cases. See the following Mathematica document for the solutions.
28
6. Kronig-Penney Model Brief note: the original work on this problem was done by Kronig and Penney in their 1931 paper, with 981 references according to Google Scholar. R. de L. Kronig and W. G. Penney. “Quantum Mechanics of Electrons in Crystal Lattices”, Proc. R. Soc. A., vol. 130, pg. 499, 1931. (a) Show that for ψ(x) = eikx u(x), the function u(x) must obey the given single particle equation. This done by replacing ψ(x) in the Schrodinger equation with the given ansatz and carrying out the partial derivatives, see the original paper. (b) Considering the given plane wave-like solutions, (i) show that they solve the differential equation for u(x), (ii) find the condition for non-trivial solutions, and (iii) that the condition can be reduced to a simpler expression given tan(θ) = tan(θ) = -
κ2 +λ2 2
(i) This is demonstrated in the paper. You can also substitute the given solutions into the equation for u(x) in both regions of the potential and verify that they are correct. (ii) The conditions for non-trivial solutions are found by obtaining the system of four equations from the boundary conditions (value and derivative continuity). Part of this involves matching the solutions u1(a-l) and u1’(a-l) to u2 (-l) and u2’(-l), which comes from the periodicity of the potential, to introduce a factor of cos(k a) into the expression. This system of equations is shown in the Kronig-Penney paper. Casting the system of equations into matrix form, Ax = 0, shows that nontrivial solutions x (i.e. nonzero coefficients) are only possible if A has a non-empty kernel. This is only true if the determinant of A is nonzero, so setting Det(A) = 0 yields the given condition. (iii) We can introduce a factor of tan(θ) by multiplying the first term in the condition by 1=cos(λl)/cos(λl), yielding the new left-hand side: LHS = Cosλ l Sinκ a -l Tan[θ]+Cosκ a -l;
This equation can be simplified somewhat using the trig. identity sin(u) tan(θ) - cos(u) = cos(u-θ)sec(θ) as follows LHS = LHS // TrigReduce// Simplify Cos[θ+(-a +l) κ] Cos[l λ] Sec[θ] κ2 +λ2
Since tan(θ) = - 2 κ λ tan(λ l) ≡ u, we can use the following identity to rewrite sec(θ)
Printed by Wolfram Mathematica Student Edition
29
2
Kronig_Penney_nontrivial.nb
Sec[ArcTan[u]] 1 +u2
This eliminates the sec(θ), yielding an expression close to what we’re aiming for. LHS = Cosκ a -l-θ Cosλ l Sqrt1 +
κ2 +λ2 2κλ
2
Tanλ l
;
The last step is to bringing the cos(λ l) term into the square root, and rewrite the resulting cos2 (θ) term 4 κ2 λ2
as cos2 (λ l) = 1 - 4 κ2 λ2 sin 2 (λ l). Slight simplification from there yields the desired condition: 2
κ2 -λ2
1 + 4 κ2 λ2 sin 2 (λ l)
1/2
cos(κ(a - l) - θ) = cos(k a)
Plots for the Kronig-Penney Model, ϵ > Vo Thanks to Kate Pistunova, who contributed these solutions. I have only made minor changes, and added some comments. The original solutions to the Kronig-Penney model can be found in the following reference. Since the paper was published in 1930, they did not have the benefit of Mathematica. They solved the problem in the limit of infinitely narrow barriers with infinite Vo, such that the product l*Vo remains constant. This is reminscent of quantum scattering, where the area (i.e. integrated potential) of the barrier is a more fundamental quantitity than either its height or its width. Taking the limit makes the problem substantially easier to study, so you have more sophisticated solutions than Kronig and Penney did! In[6]:=
k1 :=
2e
In[7]:=
l :=
In[8]:=
o := ArcTan-
In[9]:=
F1 :=
2 (e -U0) k12 +l2 2 k1 l
Tanl l1
1/2
1 +k12 -l2 2 4 k12 l2 Sinl l12
Cosk1 a -l1-o /. l1 → 0.2 /. U0 → 3 /. a → 10
This is a plot of the left-hand side of the condition for non-trivial solutions as a function of energy - since the right-hand side is cos(k a), only values in the range [-1,1] correspond to real k. In[31]:=
PlotF1 , {e, 0, 20}, ImageSize → Small 3 2
Out[31]=
1
-1
5
10
15
20
-2
In plot form, for a reasonable combination of a, l, and Vo, this results in band gaps that closely resemble those in real solids.
Printed by Wolfram Mathematica Student Edition
30
Kronig_Penney_nontrivial.nb
In[32]:=
ContourPlotCos 10 k ⩵ F1, k, 0, 1, {e, 3, 8}, GridLines→ #1 π/10, Dashed &/@Range[30], None, ImageSize → Small 8 7 6
Out[32]=
5 4 3 0.0
0.2
0.4
0.6
0.8
1.0
Since the plots are periodic in multiples of 2π/a, a given range of k values (in this case [0,1]) covers more variation when a is increased. In addition, the larger distance between barriers decreases the energy, similar to the classic infinite square well problem where the energy eigenstates decrease with increasing widths. In addition, the band gaps are smaller. This plot demonstrates the effect. In[12]:=
F12 := 1/2
2
1 +k12 -l2 4 k12 l2 Sinl l12 In[33]:=
Cosk1 a -l1-o /. l1 → 0.2 /. U0 → 3 /. a → 30
ContourPlotCos 30 k ⩵ F12, k, 0, 1, {e, 3, 8}, GridLines→ #1 π/10, Dashed &/@Range[30], None, ImageSize → Small 8 7 6
Out[33]=
5 4 3 0.0
0.2
0.4
0.6
0.8
1.0
Next we consider very small values for the width of the barriers. The band gaps decrease as shown in the following plot. In[14]:=
F13 := 2
1/2
1 +k12 -l2 4 k12 l2 Sinl l12
Cosk1 a -l1-o /. l1 → 0.05 /. U0 → 3 /. a → 10
Printed by Wolfram Mathematica Student Edition
31
3
4
Kronig_Penney_nontrivial.nb
In[34]:=
ContourPlotCos 10 k ⩵ F13, k, 0, 1, {e, 3, 8}, GridLines→ #1 π/10, Dashed &/@Range[30], None, ImageSize → Small 8 7 6
Out[34]=
5 4 3 0.0
0.2
0.4
0.6
0.8
1.0
This next plot looks at very large barrier heights Vo, where the band gaps become much larger. In[16]:=
1/2
2
F14 := 1 +k12 -l2 4 k12 l2 Sinl l12
Cosk1 a -l1-o /. l1 → 0.02 /. U0 → 300 /.
a → 10 In[35]:=
ContourPlotCos 10 k ⩵ F14, k, 0, 1, {e, 300, 350}, GridLines→ #1 π/10, Dashed &/@Range[30], None, ImageSize → Small 350 340 330
Out[35]=
320 310 300 0.0
0.2
0.4
0.6
0.8
1.0
Here are the corresponding solutions when ϵ < Vo. The plots are qualitatively similar. The same equations are used as above, but λ → ⅈλ since the argument to the square root becomes negative. This transforms the sine terms to hyperbolic sines and the cosines to hyperbolic cosines. In[18]:=
F2 := k12 -l2 Sink1 a -l1 Sinhl l1+Cosk1 a -l1 Coshl l1 /. l1 → 0.2 /. a → 10 /. U0 → 50 2 k1 l
The following is a plot of left-hand side of the new condition for non-trivial solutions. An interesting phenomenon is that the band gaps start large, nearly vanish, then become large again.
Printed by Wolfram Mathematica Student Edition
32
Kronig_Penney_nontrivial.nb
In[36]:=
Plot{F2, -1, 1}, {e, 0, 50}, PlotStyle → Automatic , DirectiveBlack, Dashed, DirectiveBlack, Dashed, PlotRange → Full, {-1.5, 1.5}, ImageSize → Small 1.5 1.0 0.5
Out[36]=
-0.5
10
20
30
40
50
-1.0 -1.5 In[37]:=
F2 := k12 -l2 Sink1 a -l1 Sinhl l1+Cosk1 a -l1 Coshl l1 /. l1 → 0.2 /. a → 10 /. U0 → 2 2 k1 l ContourPlotCos10 k ⩵ F2, k, 0, 1, {e, -0.2, 2}, GridLines→ #1 π/10, Dashed &/@Range[3], None, PlotRange → {{0, 1}, {-0.1, 2}}, ImageSize → Small 2.0
1.5
Out[38]=
1.0
0.5
0.0 0.0
0.2
0.4
0.6
0.8
1.0
In[22]:=
F22 := k12 -l2 Sink1 a -l1 Sinhl l1+Cosk1 a -l1 Coshl l1 /. l1 → 0.2 /. a → 30 /. U0 → 3 2 k1 l
In[39]:=
ContourPlotCos30 k ⩵ F22, k, 0, 1, {e, 0, 2}, GridLines→ #1 π/100, Dashed &/@Range[3], None, ImageSize → Small 2.0
1.5
Out[39]= 1.0
0.5
0.0 0.0 In[24]:=
0.2
0.4
0.6
0.8
1.0
F23 := k12 -l2 Sink1 a -l1 Sinhl l1+Cosk1 a -l1 Coshl l1 /. l1 → 0.02 /. a → 10 /. U0 → 3 2 k1 l
Printed by Wolfram Mathematica Student Edition
33
5