# Getting useful units for combinations of physical constants like on WolframAlpha?

When I need to do order of magnitude estimates, I often encounter situations where I have a mess of physical constants and difficulty putting them in sensible units.

For example, consider the combination:

$$\frac{\hbar^2}{m e}$$

Where $$\hbar$$ is the reduced Planck's constant, $$m$$ is the electron mass, and $$e$$ is the electron charge. I don't know the units of this combination a priori, but when I enter it on Wolframalpha I get really useful unit combinations (see here), and find the combination has units of $$Volt \cdot Meter^2$$.

When I put this combination into Mathematica and use

hbar = PlanckConstantReduced; m = ElectronMass; e = ElectronCharge;

UnitSimplify[(hbar ^2)/(m e)]


I get the quantity in a mess of units of Joules, Seconds, Coulombs, and Kilograms.

Is there a way to get Mathematica to suggest a list of useful units like with Wolframalpha? In other words, can we get Mathematica to give "smart" simplified unit combinations?

• I think you should use Quantity["PlanckConstantReduced"], etc. – vi pa Dec 22 '19 at 17:34
• @vipa thank you for your comment. I tried doing the same with Quantity, but it still did not give a similar output to Alpha regarding units. Instead it seemed to be trying to simplify the particular combination of constants – KF Gauss Dec 22 '19 at 17:59
• I'm sorry is Quantity["ReducedPlanckConstant"]. – vi pa Dec 22 '19 at 19:08
• UnitSimplify is supposed to do exactly what you want, but it doesn't work. Maybe file a bug report? – Roman Dec 22 '19 at 20:44

Here's an attempt at automatic unit simplification. I'll give examples first, and implementation later.

First, define some target units:

SIbase = {"Seconds", "Meters", "Kilograms", "Amperes", "Kelvins", "Moles",
SIderived = {"Newtons", "Pascals", "Joules", "Watts", "Coulombs", "Volts",
"Farads", "Ohms", "Webers", "Teslas", "Henries", "Lumens", "Lux"};
physicalconstants = {"SpeedOfLight", "PlanckConstant", "BoltzmannConstant",
"ElementaryCharge", "ElectricConstant", "MagneticConstant",
"GravitationalConstant", "JosephsonConstant"};


Express $$\hbar^2/(m e)$$ in terms of SI units:

Z = Quantity["ReducedPlanckConstant"]^2/(Quantity["ElectronMass"]*Quantity["ElectronCharge"]);
unitSimplify[Z, Join[SIbase, SIderived]]
(*    {7.61996423*10^-20 m^2 V}    *)


Express $$h/k_B^2$$ in SI units: who would have thought that the most concise way of expressing the unit is squared Kelvins per Watt?

unitSimplify[Quantity["PlanckConstant"]/Quantity["BoltzmannConstant"]^2,
Join[SIbase, SIderived]]
(*    {6626070150000000000000000/1906191661201 K^2/W}    *)


Express the speed of light in funny units: there are four equivalently simple solutions,

unitSimplify[Quantity["SpeedOfLight"],
{"Meters", "Seconds", "Hertz", "Wavenumbers"}]
(*    {299792458 m/s,
299792458 m Hz,
29979245800 per wavenumber per second,
29979245800 Hz/wavenumber}    *)


Express a meter in terms of physical constants:

unitSimplify[Quantity[1, "Meters"], physicalconstants]
(*    {2.4683*10^34 Sqrt[G] Sqrt[h]/c^(3/2)}    *)


## implementation

First, a helper function: find the sparsest vector $$\vec{x}$$ of the underdetermined linear system of equations $$\vec{x}\cdot A=\vec{b}$$: (this is an NP-complete problem used in Compressed Sensing and there are much faster heuristic algorithms available)

sparseSolve[A_, b_, n_] :=
Quiet@Check[{#, LinearSolve[Transpose[A[[#]]], b]}, Nothing] & /@
Subsets[Range[Length[A]], {n}]
sparseSolve[A_, b_] := Module[{n, solutions, shortsolutions},
(* find the solutions with smallest number n of nonzero entries *)
n = 1;
While[(solutions = sparseSolve[A, b, n]) == {}, n++];
(* of these, take the solutions that are smallest by 1-norm *)
shortsolutions = MinimalBy[solutions, Norm[#[], 1] &];
(* convert to solution vectors *)
SparseArray[Thread[Rule @@ #], Length[A]] & /@ shortsolutions]


The automated unit simplifier: simplify the units of Q by forming combinations of the units given in the list U. A list of equivalently simple solutions is returned:

unitSimplify[Q_Quantity, U_List] :=
Module[{Uunitdimensions, unitdimensionslist, Uunitexponents, Qunitexponents,
simplestunits},
(* find the unit dimensions of each unit given in the list U *)
Uunitdimensions = UnitDimensions[UnitConvert[#]] & /@ U;
(* make a list of all the unit dimensions used here *)
unitdimensionslist = Union @@ Uunitdimensions[[All, All, 1]];
(* for each unit in U, make a list of exponents of the unit dimensions in *)
(* the order of unitdimensionslist                                        *)
Uunitexponents = Lookup[Rule @@@ # & /@ Uunitdimensions, unitdimensionslist, 0];
(* for the desired unit Q, make a list of exponents of the unit dimensions *)
(* in the order of unitdimensionslist                                      *)
Qunitexponents = Lookup[Rule @@@ UnitDimensions[Q], unitdimensionslist, 0];
(* find the simplest possible ways of expressing Qunitexponents as a linear *)
(* combination of the rows of Uunitexponents                                *)
simplestunits = Times @@ (U^#) & /@ sparseSolve[Uunitexponents, Qunitexponents];
(* for each solution, convert Q to this unit *)
UnitConvert[Q, #] & /@ simplestunits]


Explicitly requesting a specific unit works fine:

hbar = Quantity["ReducedPlanckConstant"];
m = Quantity["ElectronMass"];
e = Quantity["ElectronCharge"];

UnitConvert[hbar^2/(m e), "m^2 V"]
(*    7.61996423*10^-20 m^2 V    *)


To automate the process, maybe this answer could be useful for you. Applying it to the present case, we can build a unit system that has Volts as one of the base units,

U = makeUnitSystem[{"Volts"}]
(*    {"TimeUnit" -> "Seconds", "LengthUnit" -> "Meters",
"MassUnit" -> "Kilograms", "TemperatureUnit" -> "Kelvins",
"TemperatureDifferenceUnit" -> "KelvinsDifference",
"ElectricCurrentUnit" -> ("Kilograms" ("Meters")^2)/(("Seconds")^3 "Volts"),
"LuminousIntensityUnit" -> "Candelas", "AmountUnit" -> "Moles",

(Meters are a default base unit and don't need to be mentioned). Then use this unit system U to convert the desired expression automatically to this system of base units:
unitConvert[hbar^2/(m e), U]