# Millikan's oil-drop experiment

Robert Millikan's famous oil-drop experiments were carried out at the University of Chicago from 1909 to determine the magnitude of the charge of the electron (since May 2019, this quantity has been fixed by definition at $1.602176634 \times 10^{-19}\;\mathrm{C}$.) In a single experiment, an electrically charged oil droplet was observed to fall a known distance, $d$, between two uncharged plates at its terminal velocity, $v_g$: from the time taken, $t_g$, the droplet's radius, $a$, can be deduced. Next, a voltage was applied to the plates, inducing an electric field between them. As the droplet rises under the resulting net force, the time taken, $t_e$, for it to move back up through the same distance, $d$, can be used to deduce its total charge, $q$, which is observed to be an integer multiple of the same base value, $e$, that is: $q=Ne$.

In this example we adopt a coordinate system in which the droplet's vertical position, $z$, increases in the "up" direction. For the free-fall part of the experiment, when the droplet falls at constant terminal velocity $v_g=-d/t_g$ there is no net force on it: the sum of the gravitational and drag forces is zero: $$F_g + F_d = 0 \quad \Rightarrow -m'g - 6\pi\eta a v_g = 0,$$ where $m' = \frac{4}{3}\pi a^3\rho' = \frac{4}{3}\pi a^3(\rho_\mathrm{oil} - \rho_\mathrm{air})$ is the effective mass of the droplet (after the mass of air it displaces is taken into account), $g = 9.803\;\mathrm{m\,s^{-1}}$ is the acceleration due to gravity in Chicago, and $\eta = 1.859 \times 10^{-5}\;\mathrm{kg\,m^{-1}\,s^{-1}}$ is the air viscosity under the experimental conditions (temperature, humidity, etc.). Rearranging, we get the following expression for the droplet radius: $$a = \sqrt{\frac{-9\eta v_g}{2\rho'g}}.$$ When a suitable voltage is applied to the plates and the droplet moves upwards at a constant velocity $v_e=d/t_e$, the force due to the electric field is balanced by gravity and the drag force (at this new velocity): \begin{align*} & F_e + F_g + F_\mathrm{d}' = 0 \quad \Rightarrow qE + 6\pi\eta a v_g - 6\pi\eta a v_e = 0\\ \Rightarrow \quad & q = \frac{6\pi\eta a(v_e - v_g)}{E} = \frac{6\pi\eta ad}{E}\left( \frac{1}{t_g} + \frac{1}{t_e} \right) \end{align*} Each droplet (labeled A – H) was observed three times for each different charge, $q$, acquired by exposure to X-rays (up to seven experiments per droplet).

The data below, which are also provided in the file eg10-millikan-data.txt, give the time data for a number of such experiments conducted with an oil of density $\rho_\mathrm{oil} = 917.3\;\mathrm{kg\,m^{-3}}$ on a day for which $\rho_\mathrm{air} = 1.17\;\mathrm{kg\,m^{-3}}$. The magnitude of the electric field was $E = 322.1\;\mathrm{kN\,C^{-1}}$ and the distance the drops move, $d=11.09\;\mathrm{mm}$. We can use these data to estimate $e$ (assuming it is not fixed by definition) as follows.

drop expt    tg       te       tg       te       tg       te
A       1   13.102   46.822   12.941   46.896   13.086   46.681
A       2   12.938   86.767   13.032   86.952   13.086   86.746
A       3   13.023   61.082   12.958   60.826   12.998   60.860
A       4   12.943   86.747   12.922   86.840   13.054   86.899
B       1   11.434   56.305   11.350   56.097   11.246   56.282
B       2   11.402   75.823   11.584   75.819   11.487   76.063
B       3   11.591   44.717   11.397   44.851   11.364   44.776
B       4   11.443   75.905   11.368   75.975   11.457   76.041
B       5   11.434   75.939   11.414   75.880   11.444   75.929
B       6   11.559   75.892   11.414   75.924   11.292   75.985
B       7   11.394   44.716   11.589   44.753   11.401   44.794
C       1   16.197  100.458   16.010  100.486   16.329  100.461
C       2   16.241   47.727   16.106   47.714   16.177   47.625
C       3   16.133   37.879   16.267   37.746   16.203   37.709
C       4   16.170   64.765   16.136   64.649   16.229   64.508
D       1   16.176   38.017   16.127   37.910   16.282   38.020
D       2   16.275   38.280   16.092   38.208   16.133   38.092
D       3   16.422   48.327   16.073   48.284   16.212   48.184
D       4   16.134   38.202   16.258   38.270   16.105   38.229
D       5   16.164  102.562   16.217  102.673   16.194  102.696
E       1   12.275   55.020   12.116   54.962   12.307   54.978
E       2   12.157   54.772   12.183   54.967   12.046   55.219
E       3   12.146   55.004   12.118   54.938   12.346   54.869
E       4   12.319   43.635   12.243   43.552   12.073   43.582
F       1   14.172   61.946   14.174   61.970   14.069   61.959
F       2   14.145   90.718   13.955   90.707   14.075   90.866
F       3   14.070   62.147   14.074   61.961   14.247   61.892
F       4   14.017   61.968   14.101   61.921   14.106   62.174
G       1    9.723   50.375    9.527   50.482    9.502   50.508
G       2    9.463   63.755    9.670   63.853    9.509   63.827
G       3    9.448   63.804    9.407   63.899    9.563   63.768
G       4    9.327   63.855    9.518   63.967    9.533   63.824
H       1   13.192   73.375   13.167   73.338   13.316   73.449
H       2   13.042   42.642   13.387   42.428   13.334   42.459
H       3   13.389   42.379   13.244   42.373   13.055   42.610
H       4   13.114   73.161   13.226   73.384   13.207   73.257
H       5   13.030   73.295   13.022   73.419   13.438   73.512

First, define the necessary parameters:

eta = 1.859e-5                  # air viscosity, kg.m-1.s-1
rho_air = 1.17                  # air density, kg.m-3
rho_oil = 917.3                 # oil density, kg.m-3
rhop = rho_oil - rho_air
g = 9.803                       # acceleration due to gravity, m.s-2
d = 11.09e-3                    # rise/fall distance, m
E = -322.1e3                    # electric field vector (points down!)

Next, read in the data, assigning the first two columns to a MultiIndex:

In [x]: import pandas as pd
In [x]: df = pd.read_csv('eg10-millikan-data.txt', delim_whitespace=True,
index_col=[0, 1])
Out[x]:
tg      te    tg.1    te.1    tg.2    te.2
drop expt
A    1     13.102  46.822  12.941  46.896  13.086  46.681
2     12.938  86.767  13.032  86.952  13.086  86.746
3     13.023  61.082  12.958  60.826  12.998  60.860
4     12.943  86.747  12.922  86.840  13.054  86.899
B    1     11.434  56.305  11.350  56.097  11.246  56.282

Note that pandas has added a counting integer to the column names to make them distinct.

We will start with just a single droplet, taking the transpose of its data:

In [x]: dropA = df.loc['A'].T
In [x]: dropA
Out[x]:
expt       1       2       3       4
tg    13.102  12.938  13.023  12.943
te    46.822  86.767  61.082  86.747
tg.1  12.941  13.032  12.958  12.922
te.1  46.896  86.952  60.826  86.840
tg.2  13.086  13.086  12.998  13.054
te.2  46.681  86.746  60.860  86.899

We would prefer to label each row as simply 'tg' or 'te':

In [x]: dropA.index = dropA.index.str.slice(0, 2)
In [x]: dropA
Out[x]:
expt       1       2       3       4
tg    13.102  12.938  13.023  12.943
te    46.822  86.767  61.082  86.747
tg    12.941  13.032  12.958  12.922
te    46.896  86.952  60.826  86.840
tg    13.086  13.086  12.998  13.054
te    46.681  86.746  60.860  86.899

We require the average of all of the values of $t_g$ (in the absence of the electric field the droplet takes the same time to fall the distance $d$) and the average value of $t_e$ for each column (each experiment may have a different droplet charge, but the fall – rise times are measured three times for each experiment):

In [x]: tg = dropA.loc['tg'].values.mean()
In [x]: te = dropA.loc['te'].mean()
In [x]: tg
Out[x]: 13.006916666666667
In [x]: te
Out[x]:
expt
1    46.799667
2    86.821667
3    60.922667
4    86.828667
dtype: float64

Now use the value of $t_g$ to calculate the droplet's radius:

In [x]: a = np.sqrt(9*eta*d/tg/2/rhop/g)
In [x]: a
Out[x]: 2.8181654881967875e-06

or about 2.82 μm. The charge we deduce for each experiment is:

In [x]: q = 6 * np.pi * eta * a * d / E * (1/tg + 1/te)
In [x]: q
Out[x]:
expt
1   -3.340563e-18
2   -3.005663e-18
3   -3.172143e-18
4   -3.005631e-18
dtype: float64

Repeating this for all the droplets, we can add a column, q to the DataFrame df:

for drop in df.index.levels[0]:
drop_df = df.loc[drop].T
drop_df.index = drop_df.index.str.slice(0, 2)
tg = drop_df.loc['tg'].values.mean()
te = drop_df.loc['te'].mean()
a = np.sqrt(9*eta*d/tg/2/rhop/g)
q = 6 * np.pi * eta * a * d / E * (1/tg + 1/te)
df.loc[drop, 'q'] = q.values

It is now helpful to sort the droplet charges by magnitude and to plot the sorted array and its differences:

In [x]: sorted_q = sorted(-df.loc[:, 'q'])
In [x]: plt.plot(sorted_q)
In [x]: plt.ylabel('$|q|$')
In [x]: plt.twinx()
In [x]: dq = np.diff(sorted_q)
In [x]: plt.plot(dq)
In [x]: plt.ylabel(r'$\Delta q$')
In [x]: plt.show()

It certainly seems possible that the droplet charge is always a multiple of some value between $1\times 10^{-19}\;\mathrm{C}$ and $2\times 10^{-19}\;\mathrm{C}$. We can therefore estimate the value of $|e|$:

In [x]: e_estimate = dq[(dq>1.e-19) & (dq<2.e-19)].mean()
In [x]: e_estimate
Out[x]: 1.5697150510604604e-19

We can now add a column to df for the number of elementary charges we hypothesize for each experiment:

In [x]: df['N'] = (df['q'] / e_estimate).astype(int)

Considering all the data then gives us our estimate for the magnitude of the electron charge:

In [x]: (df['q']/df['N']).mean()
Out[x]: 1.5923552150386455e-19

within 1% of the defined value.