# Definite Integration with Exponental Heaviside Step Funtions

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.

• I would introduce dimensionless quantities and consider Integrate[ r^2 (Exp[-\[Beta] v HeavisideTheta[ r - 0] HeavisideTheta[\[Sigma] - r]] - 1), {r, 0, Infinity}, Assumptions -> {\[Beta] > 0, v > 0, \[Sigma] > 0}]. You can then take the limit  v -> Infinity and recover your result. Oct 22 '19 at 15:41
• That seems like it should be a viable approach - let me tool around with it for a bit to see if I can get it to work with a symbolic approach, and I'll get back to you if I have any question. Thank you Oct 22 '19 at 15:47
• @b.gates.you.know.what I implemented this using Limit\[LetterSpace]of\[LetterSpace]Infinite\[LetterSpace]Integral = {Integrate[func_., {var_., 0, Infinity}] :> Limit[Integrate[func, {var, 0, int\[LetterSpace]var}], int\[LetterSpace]var -> Infinity]} and Simplify[Virial\[LetterSpace]B\[LetterSpace]Def /. Limit\[LetterSpace]of\[LetterSpace]Infinite\[LetterSpace]Integral /. Potential\[LetterSpace]Gradient\[LetterSpace]Def], and was able to get a good limit equation, but mathematica still is not handling the exponential heaviside step formulas in the integral. Oct 22 '19 at 16:01
• The new output takes the form of {B[Quantity[Temperature, "Kelvins"]] -> Limit[Integrate[(-1 + E^((-Infinity)*(-1 + HeavisideTheta[-Atomic\[LetterSpace]Radius + Radial\[LetterSpace]Position])* Quantity[-(1/Temperature), 1/("BoltzmannConstant"*"Kelvins")]))*Quantity[Radial\[LetterSpace]Position^2, "Picometers"^2], {Quantity[Radial\[LetterSpace]Position, "Picometers"], 0, int\[LetterSpace]var}], int\[LetterSpace]var -> Infinity]* Quantity[-2*Pi, "AvogadroConstant"]} and again produces the error mentioned above. Oct 22 '19 at 16:04
• @b.gates.you.know.what - as the above text is a bit hard to read, I have updated the picture above to include the above approach, and its output Oct 22 '19 at 16:11

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.

• This did not quite answer the question in the way that I was hoping that it would, but it got me to the solution - I had missed that the $\phi$ definition went to an indeterminate value as written and forgotten to assume that $\sigma > 0$ . Thank you for your help Oct 23 '19 at 17:51

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}]}


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