Why does the following RegionPlot command not evaluate?

Question

I have some trouble understanding the RegionPlot command. I have defined the following input

Clear["Global`*"]

volume[R_, d_] := Pi R^2 d

atomNumber[R_, d_, n_] := 2*n*volume[R, d]

\[Gamma]Fn[x_] := 
 2/x^2 (1 - Exp[-x^2] (BesselI[0, x^2] + BesselI[1, x^2]))

\[Eta]Eval[\[Lambda]_, d_, R_, rc_, n_] := 
 UnitConvert[\[Lambda]*((atomNumber[R, d, n])^2/d^2) \[Gamma]Fn[
    R/(Sqrt[2] rc)] (1 - Exp[-d^2/(4*rc^2)])]

\[Del]Fn[\[Omega]_, m_] := 
 UnitConvert[Sqrt[Quantity["ReducedPlanckConstant"]/(\[Omega]*m)]]

coherence[\[Omega]_, d_, R_, rc_, \[Lambda]_, n_, T_, m_] := 
 UnitConvert[
  4*\[Eta]Eval[\[Lambda], d, R, rc, n]*T*(\[Del]Fn[\[Omega], m])^2]

with units

Subscript[m, 0] = Quantity[1, "AtomicMassUnit"]
R = Quantity[3.6, "Micrometers"]
d = Quantity[0.25, "Millimeters"]
n = Quantity[176.2 * 10^(27), "Meters"^(-3)]
rc = Quantity[10^(-7), "Meters"]
\[Lambda] = Quantity[10^(-17), "Seconds"^(-1)]
T = Quantity[350, "Femtoseconds"]
m = Quantity[6, "AtomicMassUnit"]
\[Omega]Paper = Quantity[40, "Terahertz"]

When I evaluate my functions separately, everything works as intended:

coherence[\[Omega]Paper, d, R, rc, \[Lambda], n, T, m]
Out: 2.28009*10^-15

coherence[\[Omega]Paper, d, R, rc, \[Lambda], n, T, m] < 1
Out: True

I want to evaluate the above inequality for a range of values for $(\lambda, r_c)$. So: Define an input range for $\lambda$ and for $r_c$ $\Rightarrow$ evaluate coherence(input) < 1 $\Rightarrow$ plot all points $(\lambda, r_c)$ for which this is true.

For that, I use the following RegionPlot command:

RegionPlot[
 Evaluate[coherence[\[Omega]Paper Quantity[40, "Terahertz"], 
    d Quantity[1/4, "Millimeters"], R Quantity[18/5, "Micrometers"], 
    rc, \[Lambda], n Quantity[1762*10^(26), "Meters"^(-3)], 
    T Quantity[350, "Femtoseconds"], 
    m Quantity[6, "AtomicMassUnit"]] < 1], {rc, 
  Quantity[10^-9, "Meters"], Quantity[10^-1, "Meters"]}, {\[Lambda], 
  Quantity[10^-10, ("Seconds")^-1], Quantity[10^-1, ("Seconds")^-1]}]

However, the code doesn't even run and the notebook stops after a while. What is the issue with this code?

My thoughts:

• the numbers are too small (I can't change the numbers -> are there alternative programs to plot?)
• there is some issue with my Mathematica version strong text

Answer 1

This is a workaround that at least produces output. It produces a few error messages and took almost 1 minute to calculate, but it may be useful as a starting point. Note that the scale is log-log.

$Version  (*  12.1.0 for Linux x86 (64-bit) (March 14, 2020)  *)

RegionPlot[coherence[ωPaper, d, R,
    10^logrc  Quantity["Meters"],
    10^logλ  Quantity["Seconds"^-1],
    n, T , m ] < 1,
  {logrc, -9, -1},
  {logλ, -10, -1}] // Quiet

All error messages have been suppressed. They are underflows caused by evaluating Exp with a large negative argument. This is considered benign, since the exponential is subtracted from 1.

Another way to visualize the region by plotting the coherence as a 3D surface. An array of about 1000 data points will produce a sufficiently smooth surface. The data points can be generated and plotted like this

data = Flatten[
    Table[{logrc, log\[Lambda], coherence[\[Omega]Paper, d, R,
       10^logrc  Quantity["Meters"],
       10^log\[Lambda]  Quantity["Seconds"^-1],
       n, T , m ]},
     {logrc, -9, -1, .2},
     {log\[Lambda], -6, -1, .2}], 1]; // Quiet

ListPlot3D[data, PlotRange -> {All, All, {0, 1}}]

This plot has been cut off at z=1 in order to visualize the boundary for which coherence == 1.

The logarithmic scales in the above plots may be misleading. A linear plot, such as the one below, shows that the surface is much steeper for small values of rc than it is for the larger values. The following plot is coherence vs rc for $\lambda = 0.00002/\text{sec}$

data2 = Table[{10^logrc, coherence[\[Omega]Paper, d, R,
      10^logrc  Quantity["Meters"],
      2*^-5 /Quantity["Seconds"],
      n, T , m ]},
    {logrc, -8, -3, 0.025}]; // Quiet
ListLinePlot[data2, FrameLabel -> {"rc"}]