Definite Integration with Exponental Heaviside Step Funtions

Question

As part of an equation of state calculator that I am working on, I have to handle a virial expansion for the equation of state. As such I have to solve for the virial coefficent, $B(T)$, which is defined as

$$B(T)=-2\pi N_a\int_0^{\infty}(e^{\frac{-\phi}{\kappa_b T}}-1)R^2dR$$

Where $\phi$ is the radial potential of a given electron at a given radius away from a given atom. This function is fairly complecated, and so there are many approximations of this formula, several of which include a heavside step function. For example, if I assume that only hard sphere inteactions hold in the system

$$\phi = \infty\left(\Theta(R-0)-\Theta(R-\sigma)\right)$$

meaning that the potential is infinite for all points from radius $R=0$ to atomic radius $\sigma$, and zero for all other points. However, I am having trouble solving the integral that comes from these symbolically using Mathematica. In this case, the system should solve out as

$$\begin{align}B(T) &=-2\pi N_a\int_0^{\infty}(e^{\frac{-\infty\left(\Theta(R-0)-\Theta(R-\sigma)\right)}{\kappa_b T}}-1)R^2dR\\&= -2\pi N_a\left(\int_0^{\sigma}(e^{\frac{-\infty\left(1-0)\right)}{\kappa_b T}}-1)R^2dR + \int_{\sigma}^{\infty}(e^{\frac{-\infty\left(1-1)\right)}{\kappa_b T}}-1)R^2dR\right)\\&= -2\pi N_a\left(\int_0^{\sigma}(e^{-\infty}-1)R^2dR + \int_{\sigma}^{\infty}(e^{0}-1)R^2dR\right)\\&= -2\pi N_a\left(\int_0^{\sigma}(-1)R^2dR + \int_{\sigma}^{\infty}(0)R^2dR\right)\\&= -2\pi N_a\left(\int_0^{\sigma}(-1)R^2dR + 0\right)\\&= 2\pi N_a\int_0^{\sigma}R^2dR\\&= \frac{2\pi}{3} N_a\sigma^3 \end{align}$$

where $N_a$ is Avagadros Constant. What is important to note about this is that I have some 10 or so different sets of assumptions that I can make to define different versions of $\phi$, and that i would like to be able to solve these systems symbolically using Mathematica.

I have reduced what I have worked out so far to a near minimal example which you can see below.

$Assumptions := {Inequality[1000000, Greater, Temperature, GreaterEqual, 0], Radial\[LetterSpace]Position >= 0, Atomic\[LetterSpace]Radius >= 0, Second\[LetterSpace]Step >= 0}
Subscript[k, b] := Quantity[1, "BoltzmannConstant"]
Subscript[A, n] := Quantity[1, "AvogadroConstant"]
T := Quantity[Temperature, "Kelvins"]
R := Quantity[Radial\[LetterSpace]Position, "Picometers"]
\[Sigma] := Quantity[Atomic\[LetterSpace]Radius, "Picometers"]
Step[r_] := HeavisideTheta[QuantityMagnitude[UnitConvert[r, "Picometers"]]]
Potential\[LetterSpace]Gradient\[LetterSpace]Def = {\[Phi] -> Infinity*(Step[R + 0] - Step[R - \[Sigma]])}
Virial\[LetterSpace]B\[LetterSpace]Def = {B[T] -> -2*Pi*Subscript[A, n]*Integrate[(E^(-(\[Phi]/(Subscript[k, b]*T))) - 1)*R^2, {R, 0, Infinity}]}
Simplify[Virial\[LetterSpace]B\[LetterSpace]Def /. Potential\[LetterSpace]Gradient\[LetterSpace]Def]

This appears in the IDE as

and when executed, this produces the error

Integrate: Missing or incompatible quantities encountered in integration limits {Radial_Position pm,0,Infinity}.

and fails to produce the desired result, instead producing a function somewhere along the first line of the derivation above.

What approaches can I use to solve this type of system symbolically using Mathematica?

Approaches I have tried

I have tried restricting the bounds of several of the variables using the $Assumptions tag, declaring the relevant values as being Constants, and redefining the $\phi$ function as using piecewise notation to see if any of these help with the integration, however these approaches have as of yet proved fruitless.

I have tried to convert the integral from $0$ to $\infty$ to the limit of the integral from $0$ to $x$ as $x\to\infty$. As of yet this has also not provided an appropriate answer.

Answer 1

First of, do not use units in Mathematica code, this will only complicate things. I am sure you can keep track of units without having them spelled out within Mathematica.

Note that

$$ⅇ^{\frac{-\infty\left(\Theta(R-0)-\Theta(R-\sigma)\right)}{\kappa_b T}}=\Theta(-R)+\Theta(R-\sigma)$$

as you show in your calculation. So let us define an expression representing that

expϕ = HeavisideTheta[-r]+HeavisideTheta[r-s];

Then you can do your desired integral symbolically, first without taking the boundary values

indefiniteI = -2 π Subscript[Ν, a] Integrate[(expϕ - 1) r^2, r]

The boundary at zero is simple to evaluate:

Assuming[s > 0, indefiniteI /. r -> 0 // Simplify]

0

So we know that the result is going to be given by the boundary value at infinity. Trying to plug in r->Infinity right away confuses Mathematica, so instead we substitute in the fact that both step functions become equal to 1 as r becomes large

indefiniteI /. HeavisideTheta[_] -> 1

which is the desired result.

Alternatively, instead of directly replacing the step functions, you could expand the step function arguments to leading order in r, simplify, and then set r to infinity:

(indefiniteI /. HeavisideTheta[arg_] :> HeavisideTheta[ Series[arg, {r, Infinity, -Exponent[arg, r]}] // Normal] // Simplify)/.r->Infinity

which leads to the same result.

You could go through these steps 10 times, with updated expϕ expressions, or make a vector out of expϕ with the different cases as components and do the calculation in one go for example.

Answer 2

From the input of both b.gates.you.know, and Kagaratsch, I was able to find that the issue lied in two places with my code - namely that the definition of $\phi$ went to an indeterminite value at $R>\sigma$, and that I was failing to assume that $\sigma >0$. From this, I can get the code down to

$Assumptions := {T > 0, Subscript[k, b] > 0, \[Sigma] > 0}
Potential\[LetterSpace]Gradient\[LetterSpace]Def := {\[Phi] -> Piecewise[{{Infinity, R < \[Sigma]}, {0, R >= \[Sigma]}}]}
Virial\[LetterSpace]B\[LetterSpace]Def := {B[T] -> -2*Pi*Subscript[A, n]*Integrate[(E^(-(\[Phi]/(Subscript[k, b]*T))) - 1)*R^2, {R, 0, Infinity}]}
Virial\[LetterSpace]B\[LetterSpace]Def /. Potential\[LetterSpace]Gradient\[LetterSpace]Def

or

Potential\[LetterSpace]Gradient\[LetterSpace]Def := {\[Phi] -> Piecewise[{{Infinity, R < \[Sigma]}, {0, R >= \[Sigma]}}]}
Virial\[LetterSpace]B\[LetterSpace]Def := {B[T] -> -2*Pi*Subscript[A, n]*Integrate[(E^(-(\[Phi]/(Subscript[k, b]*T))) - 1)*R^2, {R, 0, Infinity}]}
Assuming[{R > 0, Subscript[k, b] > 0, \[Sigma] > 0}, Virial\[LetterSpace]B\[LetterSpace]Def /. Potential\[LetterSpace]Gradient\[LetterSpace]Def]

Each of which renders the desired output of

$$\left\{B[T]\rightarrow\frac{2}{3}\pi\sigma^3A_n\right\}$$

while maintaining the symbolic nature of the solution

Thank you to both @B.Gates.You.Know and @Kagaratsch for helping me find this solution