Reduced masses and frequency ratios
Clear
Clear["Global`*"]
Factors
Constants
N[I]
N[π]
N[Pi]
N[E]
[c] = m/s
c=299792458.;
pi=3.14159265; (*N[Pi,9]*)
invpi= 0.318309886; (*1/N[Pi,9]*)
Prefix Transformation
centi
base2centi=10.^(2.);
centi2base=10.^(-2.);
milli
base2milli=10.^(3.);
milli2base=10.^(-3.);
micro
base2micro=10.^(6.);
micro2base=10.^(-6.);
nano
base2nano=10.^(9.);
nano2base=10.^(-9.);
pico
base2pico=10.^(12.);
pico2base=10.^(-12.);
femto
base2femto=10.^(15.);
femto2base=10.^(-15.);
Kilo
base2kilo=10.^(-3.);
kilo2base=10.^(3.);
Mega
base2mega=10.^(-6.);
mega2base=10.^(6.);
Giga
base2giga=10.^(-9.);
giga2base=10.^(9.);
Tera
base2tera=10.^(-12.);
tera2base=10.^(12.);
Conversion factors
Angular frequency transformations
Angular frequency to frequency ([ω]= rad/s  to [ν ] = 1/s = Hz)
ω = 2πν
ν = ω/(2π)
rads2Hz= 1./(2.*pi);
Angular frequency to wavenumbers ([ω]=rad/s to [Overscript[ν, ~]]= cm^-1)
ω=2πν= 2π c Overscript[ν, ~]
Overscript[ν, ~]=ω/(2π c)
rads2wn= (1./(2.*pi*c))*(1./base2centi);
Frequency transformations
Frequency to wavenumber ([ν]=1/s=Hz  to [Overscript[ν, ~]]= cm^-1)
ν=  c Overscript[ν, ~]
Overscript[ν, ~]=ν/c
Hz2wn=(1./c)*(1./base2centi);
Frequency to angular frequency (Hz to rad/s)
ω = 2πν
Hz2rads= 2.pi;
Wavenumber transformations
Wavenumber to angular frequency (cm^-1 to rad/s)
wn2rads=2.*pi*c * (1./centi2base)
Wavenumber to frequency (cm^-1 to 1/s=Hz)
wn2Hz= c* (1./centi2base);
Formulas
Harmonic oscillator
For a harmonic oscillator, the oscillation frequency is given by ω=Sqrt[k/m], with k the force constant and m the mass. In a two particle system, the mass is replaced by the reduced mass, μ=(Subscript[m, 1]*Subscript[m, 2])/(Subscript[m, 1]+Subscript[m, 2]).
To predict a frequency change from a change in the mass of the oscillator, the ratio of the frequencies can be used: Subscript[ω, A]/Subscript[ω, B]=Sqrt[Subscript[k, A]/Subscript[μ, A]]*Sqrt[Subscript[μ, B]/Subscript[k, B]]=Sqrt[Subscript[k, A]/Subscript[k, B]]*Sqrt[Subscript[μ, B]/Subscript[μ, A]]
If the change in the force constant k is assumed to be insignificant, and as such Subscript[k, A] ≈ Subscript[k, B], this can be further simplified:
Subscript[ω, A]/Subscript[ω, B]=Sqrt[Subscript[μ, B]/Subscript[μ, A]]
redmass[m1RM_,m2RM_]:=(m1RM*m2RM)/(m1RM+m2RM);
freqratio[rm1FR_,rm2FR_]:=Sqrt[rm2FR/rm1FR];
Phonon
Optical at 0 k
Assuming a cubic crystal lattice with two different sites, the equation of motion close to the limit of k=0 is:
ω^2≈ 2C (1/Subscript[M, 1]+1/Subscript[M, 2])
The ratio for two different crystals, assuming similar force constants, is then:
Subscript[ω, A]/Subscript[ω, B]≈Sqrt[1/Subscript[M, 1,A]+1/Subscript[M, 2,A]/1/Subscript[M, 1,B]+1/Subscript[M, 2,B]]
freqratioPO0[m1aFRP_,m2aFRP_,m1bFRP_,m2bFRP_]:=Sqrt[((1/m1aFRP)+(1/m2aFRP))/((1/m1bFRP)+(1/m2bFRP))];
The acoustic phonon close to k=0 is nearly zero anyway due to a k vector dependence
Optical at π k, zone boundary
If we go to the Brillouin zone boundary for the same lattice, we find the motion to be:
ω≈Sqrt[(2C)/Subscript[M, 2]]
The ratio is then Subscript[ω, A]/Subscript[ω, B]≈Sqrt[Subscript[M, 2,B]/Subscript[M, 2,A]]if the force constants are similar. 
Subscript[M, 1] seems to be heavier usually, because then the phonon branches don't cross.
So this would be simply the root of the ratio of the anion mass in our case?
freqratioPOPI[m2aFRPP_,m2bFRPP_]:=Sqrt[m2bFRPP/m2aFRPP];
Acoustic at π k, zone boundary
If we go to the Brillouin zone boundary for the same lattice, we find the motion to be:
ω≈Sqrt[(2C)/Subscript[M, 1]]
The ratio is then Subscript[ω, A]/Subscript[ω, B]≈Sqrt[Subscript[M, 1,B]/Subscript[M, 1,A]]if the force constants are similar
We probably wouldn't see a shift in those then, because the heavier cations are extremely similar in mass.
It would be in line with the fact that our methods should excite only the optical phonons in zerost approximation due to the interaction of the motion with the electrical field. Acoustic modes would not have the charge displacement to do that.
Acoustic phonon in the middle goes to zero because of the k vector dependence.
Masses
Ion masses
mCl=35.453;
mBr=79.904;
mBu4N=242.463;
mBu3NMe=200.387;
mBu3NH=186.357;

mEt4N=130.25;
mEt3NMe=116.227;
mEt3NH=102.197;
Sqrt[mCl/mBr]
Sqrt[mBu3NMe/mBu3NH]
Sqrt[mBu4N/mBu3NMe]
Sqrt[mBu4N/mBu3NH]
Sqrt[mEt3NH/mBu3NH]
Sqrt[mEt3NMe/mBu3NMe]
Sqrt[mEt4N/mBu3NH]
Reduced masses
muBu4NCl=redmass[mCl,mBu4N]
muBu4NBr=redmass[mBr,mBu4N]
muBu3NMeCl=redmass[mCl,mBu3NMe]
muBu3NMeBr=redmass[mBr,mBu3NMe]
muBu3NHCl=redmass[mCl,mBu3NH]
muBu3NHBr=redmass[mBr,mBu3NH]

muEt4NCl=redmass[mCl,mEt4N]
muEt3NMeCl=redmass[mCl,mEt3NMe]
muEt3NHCl=redmass[mCl,mEt3NH]
Frequency Ratios Harmonic oscillator
Anion Change Cl-Br
frBu4N=freqratio[muBu4NCl,muBu4NBr]
frBu3NMe=freqratio[muBu3NMeCl,muBu3NMeBr]
frBu3NH=freqratio[muBu3NHCl,muBu3NHBr]
1/frBu3NMe
1/frBu3NH
1/frBu4N
Cation Change Cl
frClBu4NBu3NMe=freqratio[muBu4NCl,muBu3NMeCl]
frClBu4NBu3NH=freqratio[muBu4NCl,muBu3NHCl]
frClBu3NMeBu3NH=freqratio[muBu3NMeCl,muBu3NHCl]

frClBu4NEt4N=freqratio[muBu4NCl,muEt4NCl]
frClBu3NMeEt3NMe=freqratio[muBu3NMeCl,muEt3NMeCl]
frClBu3NHEt3NH=freqratio[muBu3NHCl,muEt3NHCl]
Actually took the inverses (divided by Bu4N frequency) bc I took the one that gave a ratio above 1, so divided by lower frequency.
frClBu3NHBu4N=1/frClBu4NBu3NH
frClBu3NMeBu4N=1/frClBu4NBu3NMe
frClBu3NHBu3NMe=1/frClBu3NMeBu3NH

frClEt4NBu4N= 1/frClBu4NEt4N
frClEt3NMeBu3NMe=1/frClBu3NMeEt3NMe
frClEt3NHBu3NH=1/frClBu3NHEt3NH
Cation Change Br
frBrBu4NBu3NMe=freqratio[muBu4NBr,muBu3NMeBr]
frBrBu4NBu3NH=freqratio[muBu4NBr,muBu3NHBr]
frBrBu3NMeBu3NH=freqratio[muBu3NMeBr,muBu3NHBr]
Actually took the inverses (divided by Bu4N frequency) bc I took the one that gave a ratio above 1, so divided by lower frequency.
frBrBu3NHBu4N=1/frBrBu4NBu3NH
frBrBu3NMeBu4N=1/frBrBu4NBu3NMe
frBrBu3NHBu3NMe=1/frBrBu3NMeBu3NH
Frequency Ratios Phonons Zone boundary
Anion Change Cl-Br
Only one because cation mass change doesn't influence the ratio
frzb=freqratioPOPI[mCl,mBr]
1/frzb
Cation Change Cl
Shouldn' t see the cation change bc anion, which should be the only determining factor, stays the same
Cation Change Br
Shouldn' t see the cation change bc anion, which should be the only determining factor, stays the same
Frequency Ratios Phonons middle/0k
Anion Change Cl-Br
fr0kBu4N=freqratioPO0[mBu4N,mCl,mBu4N,mBr]
fr0kBu3NMe=freqratioPO0[mBu3NMe,mCl,mBu3NMe,mBr]
fr0kBu3NH=freqratioPO0[mBu3NH,mCl,mBu3NH,mBr]
1/fr0kBu3NH
1/fr0kBu4N
1/fr0kBu3NMe
Cation Change Cl
fr0kBu3NHBu4N=freqratioPO0[mBu3NH,mCl,mBu4N,mCl]
fr0kBu3NMeBu4N=freqratioPO0[mBu3NMe,mCl,mBu4N,mCl]
fr0kBu3NHBu3NMe=freqratioPO0[mBu3NH,mCl,mBu3NMe,mCl]

fr0kEt4NBu4N=freqratioPO0[mEt4N,mCl,mBu4N,mCl]
fr0kEt3NMeBu3NMe=freqratioPO0[mEt3NMe,mCl,mBu3NMe,mCl]
fr0kEt3NHBu3NMH=freqratioPO0[mEt3NH,mCl,mBu3NH,mCl]
1/fr0kBu3NHBu4N
1/fr0kBu3NMeBu4N
1/fr0kBu3NHBu3NMe

1/fr0kEt4NBu4N
1/fr0kEt3NMeBu3NMe
1/fr0kEt3NHBu3NMH
Cation Change Br
fr0kBu3NHBu4N=freqratioPO0[mBu3NH,mBr,mBu4N,mBr]
fr0kBu3NMeBu4N=freqratioPO0[mBu3NMe,mBr,mBu4N,mBr]
fr0kBu3NHBu3NMe=freqratioPO0[mBu3NH,mBr,mBu3NMe,mBr]
Force Constants
The force constant of a motion (if approximated as a harmonic oscillator) can be calculated as k=mω^2
k1OVERk2[m1_,o1_,m2_,o2_]:=(m1* o1^2)/(m2* o2^2);
Compare this to the intensity changes expected from the mass change after Subscript[μ, 01]^2=planck/2 1/(√m*k) (δμ/δq)^2=A which would be Subscript[A, Br]/Subscript[A, Cl]=√((Subscript[m, Cl] Subscript[k, cl])/(Subscript[m, Br] Subscript[k, Br]))
Frequencies
fABu3NHBr=38.46;
fCBu3NHBr=121.18;
fABu3NHCl=61.21;
fCBu3NHCl=174.78;

fABu3NMeBr=25.49;
fCBu3NMeBr=77.62;
fABu3NMeCl=45.14;
fCBu3NMeCl=105.78;

fABu4NBr=32.53;
fCBu4NBr=84.47;
fABu4NCl=47.13;
fCBu4NCl=88.99;
Amplitudes
AABu3NHBr=0.97;
ACBu3NHBr=0.56;
AABu3NHCl=0.77;
ACBu3NHCl=0.62;

AABu3NMeBr=0.73;
ACBu3NMeBr=0.33;
AABu3NMeCl=0.66;
ACBu3NMeCl=0.53;

AABu4NBr=17.52;
ACBu4NBr=0.46;
AABu4NCl=7.26;
ACBu4NCl=2.87;
kCl/kBr - Reduced Mass
Bu3NH
kArmBu3NHCl=k1OVERk2[muBu3NHCl,fABu3NHCl,muBu3NHBr,fABu3NHBr]
kCrmBu3NHCl=k1OVERk2[muBu3NHCl,fCBu3NHCl,muBu3NHBr,fCBu3NHBr]
Bu3NMe
kArmBu3NMeCl=k1OVERk2[muBu3NMeCl,fABu3NMeCl,muBu3NMeBr,fABu3NMeBr]
kCrmBu3NMeCl=k1OVERk2[muBu3NMeCl,fCBu3NMeCl,muBu3NMeBr,fCBu3NMeBr]
Bu4N
kArmBu4NCl=k1OVERk2[muBu4NCl,fABu4NCl,muBu4NBr,fABu4NBr]
kCrmBu4NCl=k1OVERk2[muBu4NCl,fCBu4NCl,muBu4NBr,fCBu4NBr]
kCl/kBr - Anion Mass
Bu3NH
k1OVERk2[mCl,fABu3NHCl,mBr,fABu3NHBr]
k1OVERk2[mCl,fCBu3NHCl,mBr,fCBu3NHBr]
Bu3NMe
k1OVERk2[mCl,fABu3NMeCl,mBr,fABu3NMeBr]
k1OVERk2[mCl,fCBu3NMeCl,mBr,fCBu3NMeBr]
Bu4N
k1OVERk2[mCl,fABu4NCl,mBr,fABu4NBr]
k1OVERk2[mCl,fCBu4NCl,mBr,fCBu4NBr]
ABr/ACl 
Bu3NH
measured
AABu3NHBr/AABu3NHCl
ACBu3NHBr/ACBu3NHCl
theoretical - reduced mass
Sqrt[(kArmBu3NHCl*(muBu3NHCl/muBu3NHBr))]
Sqrt[(kCrmBu3NHCl*(muBu3NHCl/muBu3NHBr))]
(kArmBu3NHCl*(muBu3NHCl/muBu3NHBr))
(kCrmBu3NHCl*(muBu3NHCl/muBu3NHBr))
theoretical - anion mass
Sqrt[(kArmBu3NHCl*(mCl/mBr))]
Sqrt[(kCrmBu3NHCl*(mCl/mBr))]
(kArmBu3NHCl*(mCl/mBr))
(kCrmBu3NHCl*(mCl/mBr))
Bu3NMe
measured
AABu3NMeBr/AABu3NMeCl
ACBu3NMeBr/ACBu3NMeCl
theoretical - reduced mass
Sqrt[(kArmBu3NMeCl*(muBu3NMeCl/muBu3NMeBr))]
Sqrt[(kCrmBu3NMeCl*(muBu3NMeCl/muBu3NMeBr))]
(kArmBu3NMeCl*(muBu3NMeCl/muBu3NMeBr))
(kCrmBu3NMeCl*(muBu3NMeCl/muBu3NMeBr))
theoretical - anion mass
Sqrt[(kArmBu3NMeCl*(mCl/mBr))]
Sqrt[(kCrmBu3NMeCl*(mCl/mBr))]
(kArmBu3NMeCl*(mCl/mBr))
(kCrmBu3NMeCl*(mCl/mBr))
Bu4N
measured
AABu4NBr/AABu4NCl
ACBu4NBr/ACBu4NCl
theoretical - reduced mass
Sqrt[(kArmBu4NCl*(muBu4NCl/muBu4NBr))]
Sqrt[(kCrmBu4NCl*(muBu4NCl/muBu4NBr))]
(kArmBu4NCl*(muBu4NCl/muBu4NBr))
(kCrmBu4NCl*(muBu4NCl/muBu4NBr))
theoretical - anion mass
Sqrt[(kArmBu4NCl*(mCl/mBr))]
Sqrt[(kCrmBu4NCl*(mCl/mBr))]
(kArmBu4NCl*(mCl/mBr))
(kCrmBu4NCl*(mCl/mBr))
Comparison with Ludwig Paper
Additional masses
mCH3SO3=94.98;
mCF3SO3=148.95;
mNTf2=279.92;
mC2mim=111.09;
Frequencies
fEt3NHCH3SO3=149.4;
fEt3NHCF3SO3=128.9;
fEt3NHNTf2=105.2;
fC2mimNTf2=83.2;

fBu3NHCl=174.78;
fBu3NHCBr=121.18;

fBu4NCl=88.99;
fBu4NBr=84.47;
Shifts zone center
Assuming a cubic crystal lattice with two different sites, the equation of motion close to the limit of k=0 is:
ω^2≈ 2C (1/Subscript[M, 1]+1/Subscript[M, 2])
The ratio for two different crystals, assuming similar force constants, is then:
Subscript[ω, A]/Subscript[ω, B]≈Sqrt[1/Subscript[M, 1,A]+1/Subscript[M, 2,A]/1/Subscript[M, 1,B]+1/Subscript[M, 2,B]]
freqratioPO0[m1aFRP_, m2aFRP_, m1bFRP_, m2bFRP_] := Sqrt[((1/m1aFRP) + (1/m2aFRP))/((1/m1bFRP) + (1/m2bFRP))];
Calculated
freqratioPO0[mCH3SO3, mEt3NH,mCl, mBu3NH]
freqratioPO0[mCF3SO3, mEt3NH,mCl, mBu3NH]
freqratioPO0[mNTf2, mEt3NH,mCl, mBu3NH]
freqratioPO0[mCH3SO3, mEt3NH,mBr, mBu3NH]
freqratioPO0[mCF3SO3, mEt3NH,mBr, mBu3NH]
freqratioPO0[mNTf2, mEt3NH,mBr, mBu3NH]
freqratioPO0[mNTf2, mC2mim,mCl, mBu4N]
freqratioPO0[mNTf2, mC2mim,mBr, mBu4N]
freqratioPO0[mNTf2, mEt3NH,mCH3SO3, mEt3NH]
freqratioPO0[mNTf2, mC2mim,mCH3SO3, mEt3NH]
freqratioPO0[mNTf2, mEt3NH,mNTf2, mC2mim]
Experimental
Divide by lighter anion frequency
fEt3NHCH3SO3/fBu3NHCl
fEt3NHCF3SO3/fBu3NHCl
fEt3NHNTf2/fBu3NHCl
fEt3NHCH3SO3/fBu3NHCBr
fEt3NHCF3SO3/fBu3NHCBr
fEt3NHNTf2/fBu3NHCBr
fC2mimNTf2/fBu4NCl
fC2mimNTf2/fBu4NBr
fEt3NHNTf2/fEt3NHCH3SO3
fC2mimNTf2/fEt3NHCH3SO3
fEt3NHNTf2/fC2mimNTf2
Shifts zone edge
If we go to the brillouin zone boundary for the same lattice, we find the motion to be:
ω≈Sqrt[(2C)/Subscript[M, 2]]
The ratio is then Subscript[ω, A]/Subscript[ω, B]≈Sqrt[Subscript[M, 2,B]/Subscript[M, 2,A]]if the force constants are similar. 
Not quite sure how to assign Subscript[M, 1] and Subscript[M, 2,] though Subscript[M, 1] seems to be heavier usually, because then the phonon branches don't cross.
So this would be simply the root of the ratio of the anion mass in our case?
freqratioPOPI[m2aFRPP_, m2bFRPP_] := Sqrt[m2bFRPP/m2aFRPP];
Calculated
freqratioPOPI[mCH3SO3, mCl]
freqratioPOPI[mCF3SO3, mCl]
freqratioPOPI[mNTf2, mCl]
freqratioPOPI[mCH3SO3, mBr]
freqratioPOPI[mCF3SO3, mBr]
freqratioPOPI[mNTf2, mBr]
freqratioPOPI[mNTf2, mCH3SO3]
Experimental
Divide by lighter anion frequency
fEt3NHCH3SO3/fBu3NHCl
fEt3NHCF3SO3/fBu3NHCl
fEt3NHNTf2/fBu3NHCl
fEt3NHCH3SO3/fBu3NHCBr
fEt3NHCF3SO3/fBu3NHCBr
fEt3NHNTf2/fBu3NHCBr
fC2mimNTf2/fBu4NCl
fC2mimNTf2/fBu4NBr
fEt3NHNTf2/fEt3NHCH3SO3
fC2mimNTf2/fEt3NHCH3SO3
fEt3NHNTf2/fC2mimNTf2