Assumptions on unknown functions

Question

I'm trying to do analytic calculations in a quantum mechanic harmonic oscillator basis. Specifically I want to be able to evaluate functions of the many particle density.

I define the following functions: [definitions]

(* Harmonic Oscillator Eigenstate *)
basis[i_, z_] = Pi^(-1/4)/Sqrt[2^i Factorial[i]] Exp[-z^2/2] HermiteH[i, z];
(* Orbital *)
orb[i_, z_] = Sum[c[i, k] basis[k, z], {k, 0, t}];
(* Particle density of one orbital *)
orbitaldens[i_, z_] = Conjugate[orb[i, z]] psi[i, z];
(* Total density *)
dens[z_] = Sum[a[i] orbitaldens[i, z], {i, 1, n}];

To help Mathematica I define some assumptions beforehand: [assumptions]

$Assumptions = Element[z, Reals] && (* Real-space coordinate is real *)
               Element[{n, t}, Integers] && n > 0 && t > 0 && (* Just indices *)
               Element[c[i_, j_], Complexes] && (* Complex state vector *)
               Element[a[i_], Reals]; (* Real filling factor *)

What is the correct way to define assumptions on an unknown function? Here c[i_, j_] should take integers i, and j and return a complex number. a[i_] on the other hand should take on integer i and return a positive real number.

Unfortunately, the evaluation of the definitons takes very long if I include these assumptions. Without them the definitions block takes about a second to evaluate.

However, I need those assumptions, so Mathematica can do useful simplifications. For example Conjugate[orb[i,z]] orb[i,z]//FullSimplify should simplify to Abs[orb[i,z]]^2. But the same should work, if I only apply the complex conjugation to c[i,k] inside orb[i,z], because basis[i,z] is real. So far, this takes too long to evaluate.

Is it possible to speed this whole thing up?

Another problem I found is the codomain of the Factorial The following two expressions do not evaluate to True, neither to False. They just repeat the first parameter of Refine:

Refine[Element[Sqrt[Factorial[i]], Reals], Element[i, Integers] && i > 0]
Refine[Factorial[i] > 0, Element[i, Integers] && i > 0]

In words, the factorial of a positive integer is not necessarily positive.

What is wrong there? I don't see how the factorial of a positive integer could become negative. (In symbolic maths)

Answer 1

A somewhat related question may be worth linking:

Attaching persistent assumptions to symbol definition

Maybe the following will help with the problem of how to specify assumptions about coefficients in an expansion. I did something like this in another answer, and it could be applied here to your coefficients a[i]. The idea is that once the dimension of the arrays, n, is fixed, we can make dummy arrays of that specific length for the purposes of further symbolic manipulations, and impose the reality "assumption" by using an UpValue definition.

The naming of the individual array elements works by using Unique which automatically picks an unused name (starting with "a" in this case).

n = 3;

a = Table[Unique["a"], {n}];
Map[(# /: Conjugate[#] = #) &, a];

The last line is the crucial definition, saying simply that each element of the array a is equal to its conjugate (hence: real). Now I define the superposition function with x as the variable:

f = a.Sin[2 Pi Range[n] x]

(* ==> a29 Sin[2 Pi x] + a30 Sin[4 Pi x] + a31 Sin[6 Pi x] *)

For simplicity, I chose trigonometric functions as the basis.

Before doing the absolute square, let's define another almost identical function for comparison:

b = Table[Unique["b"], {n}];

g = b.Sin[2 Pi Range[n] x]

(* ==> b32 Sin[2 Pi x] + b33 Sin[4 Pi x] + b34 Sin[6 Pi x] *)

Here I left out the reality assumption for b.

For the last function, the modulus for real argument x is

Assuming[{x ∈ Reals}, Simplify[Conjugate[g] g]]

(*
==> Conjugate[
  b32 Sin[2 Pi x] + b33 Sin[4 Pi x] + 
   b34 Sin[6 Pi x]] (b32 Sin[2 Pi x] + b33 Sin[4 Pi x] + 
   b34 Sin[6 Pi x])
*)

The Conjugate remains in the result. On the other hand, if we do the same for the first function, the reality of a is used:

Assuming[{x ∈ Reals}, Simplify[Conjugate[f] f]]

(*
==> (a29 Sin[2 Pi x] + a30 Sin[4 Pi x] + 
  a31 Sin[6 Pi x])^2
*)

Edit

The reason I associated the realness condition with a using TagSet instead of using assumptions is that some tests don't recognize assumptions: see for example

Assuming[{c ∈ Reals}, Simplify[HermitianMatrixQ[{{c}}]]]

which should yield True but doesn't (because it does a literal comparison). This may not be a problem in your application, but it can't hurt to try and avoid such issues by associating the desired properties more closely with the symbol representing your coefficients or functions.

As an alternative, you could make the above definition work with the coefficients b by defining the function

h = Re[b].Sin[2 Pi Range[n] x]

where Re is put in explicitly to wrap the coefficients.