LOG#148. Path integral (III).

Posted on by December 6, 2019

Liekato

Round 3! Fight… with path integrals! XD

Introduction: generic aspects

Consider a particle moving in one dimension (the extension to ND is trivial), the hamiltonian being of the usual form:

H=\dfrac{p^2}{2m}+V(q)

The fundamental question and problem in the path integral (PI) approach of quantum mechanics (QM) is:

If the particle is at a position q (fluctuating, of course) at a time t=0, what is the probability amplitude that it will be at some other position, say q', at a time t=T?

Any other “amplitude” can be written in terms of the previous and fundamental amplitude above. It is easy to get a formal expression for this amplitude in terms of the usual formulation of QM simply solving the Schrödinger equation. Furthermore, you can make a guess using the evolution operator formalism as well! Dirac did it himself, but we had to wait for Feynman to “clean” and understand the formalism. Let me introduce the eigenstates of the position operator \hat{q}. It forms what physicists and mathematicians call “complete, orthonormal set”. Without entering into the technical (and axiomatic) issues, we have:

\hat{q}\vert q\rangle=q\vert q\rangle

\langle q'\vert q\rangle=\delta (q'-q)

\int dq \vert q \rangle \langle q\vert=1

Then, the initial state is

\vert \Psi (0)\rangle=\vert q\rangle

Letting the state evolve in time and projecting on the state \vert q'\rangle, we have the amplitude

    \[ \boxed{A(q',q;T)=\langle q'\vert \Psi (T)\rangle\equiv K(q',q;T,0)=\langle q'\vert e^{-iHT}\vert q\rangle}\]

when there is the possibility of an ambiguity, operators will be denoted with a “hat”. Otherwise, we will drop that symbol. We also fix the rationalized Planck constant \hbar to the unit, i.e., we write \hbar=1.

The object K(q',q;T,0) is called the propagator (sometimes Green function) from the initial spacetime point (q,0) to the final point (q',T). The propagator is clearly independent of the origin of time, for time independent hamiltonians:

K(q',q;T+t,t)=K(q',q;t,0)

We will derive an expression for this amplitude in the form of a sum over all possible path between the initial and final points. For this purpose, let us separate the time evolution in the above amplitude into two smaller time evolutions, with the aid of the exponential function:

e^{-iHT}=e^{-iH(T-t_1)}e^{-iHt_1}

Thus, the amplitude becomes

A=\langle q'\vert e^{-iH(T-t_1)}e^{-iHt_1}\vert q\rangle

Inserting a factor “1” in the form of the previously seen integrals for an orthonormal set, we get the next sum over the position eigenstates:

\displaystyle{A=\langle q'\vert\int dq_1\vert q_1\rangle \langle q_1\vert e^{-iHt_1}\vert q\rangle=\int dq_1 K(q',q_1;T,t_1)K(q_1,q;t_1,0)}

This formula is just an expression of the quantum mechanical rule for combining amplitudes: if a process can occur a number of ways, the amplitudes for each of these will add. And since amplitude is a complex number (!), therefore interference (constructive and destructive) will arise at some points in general. In particular, any particle, in its propagation from q to q', must be somewhere at an intermediate time t_1. However, this “somewhere” is a bit “fuzzy”, since it is spread along every possible QM path, not only the classical trajectories are possible. Of course, the weights of paths are not equal in general, but, in principle, every path does contribute! Labeling an intermediate position with q_1, we can compute the amplitude for the propagation via the point q_1, and it is the product of the two propagators as it is shown above. In general, we have always integrated over all possible intermediate (generally non-observable and virtual) positions. This seems a bit of mathematical trickery but it is all true. Quantum particles “smell” every possible path. We can repeat the division of the time interval T; let me divide it up into a large number, say N, of time intervals of tiny duration \delta=T/N. Then, we can write for the propagator the next formula:

A=\langle q'\vert (e^{-iH\delta})^{N}\vert q\rangle=\langle q'\vert e^{-iH\delta}\cdots (N-times)\cdots e^{-iH\delta}\vert q\rangle

We can again insert some complete sets of states between every exponential function, and it yields

\displaystyle{A=\langle q'\vert e^{-iH\delta}\int dq_{N-1}\vert q_{N-1}\rangle\langle q_{N-1}\vert e^{-iH\delta}\int dq_{N-2}\rangle q_{N-2}\langle q_{N-2}\vert e^{-iH\delta}\cdots \int dq_1\vert q_1\rangle\langle q_1\vert e^{-iH\delta}\vert q\rangle}

We can rewrite this expression as follows:

\displaystyle{A=\int dq_1\cdots dq_{N-1}\langle q'\vert e^{-iH\delta}\vert q_{N-1}\rangle\langle q_{N-1}\vert e^{iH\delta}\vert q_{N-2}\rangle\cdots\langle q_1\vert e^{-iH\delta}\vert q\rangle}

or, equivalently,

\displaystyle{\boxed{A=\int dq_1\cdots dq_{N-1} K_{q_{N},q_{N-1}}K_{q_{N-1},q_{N-2}}\cdots K_{q_2,q_1}K_{q_1,q_0}=\int Dq K(q_N,q_0)}}

with

K_{q_N,q_0}=K_{q_{N},q_{N-1}}\cdots K_{q_2,q_1}K_{q_1,q_0}

and

\displaystyle{\int Dq=\int d^{N-1}q=\int dq_1\cdots dq_{N-1}}

and where we have defined q_0=q and q'=q_N. Note that these initial and final positions are NOT integrated over, but only the intermediate states. This expression says that the amplitude IS the integral of the amplitude of all N-leg paths. Apart from some additional mathematical tools and details concerning the limit when N\rightarrow \infty, this is clearly going to become a sum over every possible path of the amplitude for each path. Namely, we have

A=\displaystyle{\sum_{paths} A_{path}}

where at formal level, we identify

\displaystyle{\sum_{paths}=\int Dq=\int dq_1\cdots dq_{N-1}}

and

\displaystyle{A_{path}=K_{q_N,q_0}=K_{q_{N},q_{N-1}}\cdots K_{q_2,q_1}K_{q_1,q_0}}

A final and important result is the so-called Lie-Kato-Trotter formula. For any pair of matrices/operators A,B, under very general conditions and assumptions, we have the identity:

    \[ \displaystyle{\boxed{e^{A+B}=\lim_{N\rightarrow \infty} \left(e^{A/N}e^{B/N}\right)^N}}\]

Let us make use of the Lie-Kato-Trotter formula to derive the so-called Dirac formula. For the hamiltonian

H=H_0+V(q)

with H_0=\dfrac{p^2}{2} (mass equal to one) and \delta=\dfrac{t}{N}

\displaystyle{e^{-iHt}=\lim_{N\rightarrow\infty}\left(e^{-iH_0\delta}e^{-iV\delta}\right)^N=\lim_{N\rightarrow\infty} U_\delta^N}

where

U_\delta=e^{-iH_0\delta}e^{-iV\delta}

The amplitude reads

A_\delta=\langle q'\vert U_\delta\vert q\rangle=\langle q'\vert e^{-iH_0\delta}\vert q\rangle e^{-iV(q)\delta}

We insert the set

\displaystyle{1=\int dp\vert p\rangle \langle p\vert}

and thus

\displaystyle{A_\delta=\langle q'\vert e^{-iH_0\delta}\int dp\vert p\rangle\langle p\vert \vert q\rangle e^{-iV(q)\delta}}

Carrying out the operations, we get

\displaystyle{A_\delta=\int dp \langle q'\vert p\rangle \langle p\vert q\rangle e^{-ip^2/2\delta}e^{-iV(q)\delta}}

or

\displaystyle{A_\delta=\int dp\langle q'\vert p\rangle\langle p\vert q\rangle e^{-i\left(\dfrac{p^2}{2}+V(q)\right)\delta}}

And now, using

\langle q'\vert p\rangle \langle p\vert q\rangle=\dfrac{1}{2\pi}e^{ip(q'-q)\delta}

\displaystyle{\int dp e^{ipx}e^{ip^2\delta/2}=e^{-x^2/2\delta}}

we finally obtain

\displaystyle{A_\delta=\int dp e^{ip(q'-q)}e^{-iH(\frac{p^2}{2}+V(q))\delta}}

namely,

A_\delta=\dfrac{1}{\sqrt{2\pi i \delta}}e^{-i(q'-q)^2/2\delta-iV(q)\delta}=\dfrac{1}{\sqrt{2\pi i\delta}}e^{i\delta\left(\frac{1}{2}\left(\frac{q'-q}{\delta}\right)^2-V(q)\right)}

i.e., the modern language version of the Dirac formula in terms of the classical action reads as follows

\boxed{A_\delta=\dfrac{1}{\sqrt{2\pi i\delta}}e^{i\delta S}}

Now we are ready to writhe the propagator and the PI in the configuration space. Let me write

\displaystyle{A=\langle q'\vert e^{-iHt}\vert q\rangle=\lim_{N\rightarrow\infty} \langle q'=x_N\vert U_\delta^N\vert q\rangle=x_0}

or

A=\displaystyle{\lim_{N\rightarrow\infty}\int dx_{N-1}\cdots dx_1\langle q'\vert U_\delta\vert x_{N-1}\rangle\langle x_{N-1}\vert U_\delta \vert x_{N-2}\rangle\cdots\langle x_1\vert U_\delta\vert q\rangle}

We apply now the above Dirac formula:

\displaystyle{A=\lim_{N\rightarrow\infty}\dfrac{1}{C_\delta^N}\int dx_{N-1}\cdots dx_1e^{i\Delta S(q',x_{N-1})}e^{i\Delta S(x_{N-1},x_{N-2})}\cdots e^{i\Delta S(x_1,x_0)}}

\displaystyle{\boxed{A=\lim_{N\rightarrow\infty}\dfrac{1}{C_\delta^N}\int dx_{N-1}\cdots dx_1 e^{i\Delta S(q',x_{N-1})+\cdots+i\Delta S(x_1,x_0)}=\lim_{N\rightarrow \infty}\left(\dfrac{2\pi i t}{N}\right)^{N/2}\int dx_{N-1}\cdots dx_1e^{iS(q',x_{N-1},\ldots,x_1,x_0)}}}

and where

\displaystyle{S(q',x_{N-1},\ldots,x_1,x_0)=\sum_{j}^{N-1}\dfrac{t}{N}\left(\left(\dfrac{q_{j+1}-q_j}{t/N}\right)^2-V(q_j)\right)}

The argument inside of the exponential function is a discrete approximation of the action of a path passing through the points q_0=q,q_1,q_2,\ldots,q_N=q'. Moreover, in the continuum limit, we obtain the formal definition of a Feynman PI as follows:

\displaystyle{\boxed{A=\langle q'\vert U(t,t_0)\vert q\rangle=\int \mathcal{D}qe^{iS(q)}}}

and where q_0=x_0=q(t_0) and q=q(t)=x are the initial and the final points in the amplitude. This final result is the configuration space path integral. It should be viewed as a notation of an infinite dimensional integral, or infinite differential form, of the more precise discrete expression as N\rightarrow\infty. Furthermore, this formula says that any particle, in going from one position to another, takes ALL possible path between these two points. It can be shown that, while all paths contribute, the classical path is the dominant one and it minimizes (or maximizes) the path integral. A key observation for these facts is to consider how different and close paths interfere to each other.

We can perform all this formalism in the phase space as well! Take

K_{q_{j+1},q_j}=\langle q_{j+1}\vert e^{-iH\delta}\vert q_j\rangle

as the propagator for one subinterval. We can expand the exponential function for small \delta

K_{q_{j+1},q_j}=\langle q_{j+1}\vert \left(1-iH\delta-\dfrac{1}{2}H^2\delta^2+\cdots\right)\vert q_j\rangle

It yields

K_{q_{j+1},q_j}=\langle q_{j+1}\vert q_j\rangle-i\delta\langle q_{j+1}\vert H\vert q_j\rangle+o(\delta^2)

The first term is a delta function, and we write

\langle q_{j+1}\vert q_j\rangle=\delta (q_{j+1}-q_j)=\displaystyle{\int\dfrac{dp}{2\pi}e^{ip_j(q_{j+1}-q_j)}}

The second term can be worked out as well inserting a complete set of momentum eigenstates between H and \vert q_j\rangle. This gives:

-i\delta \langle q_{j+1}\vert \left(\dfrac{p^2}{2m}+V(q)\right)\int \dfrac{dp_j}{2\pi}\vert p_j\rangle \langle p_j\vert q_j\rangle=-i\delta\int \dfrac{dp_j}{2\pi}\left(\dfrac{p^2_j}{2m}+V(q_{j+1})\right)\langle q_{j+1}\vert p_j\rangle\langle p_j\vert q_j\rangle

or equivalently

\displaystyle{-i\delta \int\dfrac{dp_j}{2\pi}\left(\dfrac{p_j^2}{2m}+V(q_{j+1})\right)e^{ip_j(q_{j+1}-q_j)}}

where we used that \langle q\vert p\rangle=\exp (ipq). We must view the momentum in the hamiltonian as the momentum operator operating to the right, while the potential (operator valued!) should operate to the left. Our last expression is asymmetric between q_j and q_{j+1}, but the origin of this is our choice of putting the factor “1” to the right of the hamiltonian H. Had we write it to the left instead to the right of the hamiltonian, we would obtain the potential V(q_j) instead V(q_{j+1}). If we introduce the shift

\overline{q_j}=\dfrac{1}{2}\left(q_j+q_{j+1}\right)

and thus V(\overline{q_j}), the expression would be symmetric in this variable. The exact choice does NOT matter in the continuum limit of a path integral! Now, we calculate

\displaystyle{K_{q_{j+1},q_j}=\int\dfrac{dp_j}{2\pi}e^{ip_j(q_{j+1}-q_j)}\left(1-i\delta\left(\dfrac{p_j^2}{2m}+V(\overline{q_j})\right)+o(\delta^2)\right)}

and

l \displaystyle{K_{q_{j+1},q_j}=\int\dfrac{dp_j}{2\pi}e^{ip_j(q_{j+1}-q_j)}e^{-i\delta H(p_j,\overline{q_j})}(1+o(\delta^2))}

There, we have N such factors in the amplitude. We can combine them, and write

\dot{q}=\dfrac{dq}{dt}=\dfrac{q_{j+1}-q_j}{\delta}

and thus

\displaystyle{A_{path}=\int\prod_{j=0}^{N-1}\dfrac{dp_j}{2\pi}e^{i\delta}\sum_{j=0}^{N-1}\left(p_j\dot{q_j}-H(p_j,\overline{q_j})\right)}

and where we have neglected a multiplicative factor

(1+o(\delta^2))^N

which will tend to one in the continuum limit N\rightarrow\infty. Then, the propagator finally becomes

\displaystyle{\boxed{K=\int dq_1\cdots dq_{N-1}A_{path}=\int\prod_{j=1}^{N-1}dq_j\int\prod_{j=0}^{N-1}\dfrac{dp_j}{2\pi}e^{i\delta}\sum_{j=0}^{N-1}\left(p_j\dot{q_j}-H(p_j,\overline{q_j}\right)}}

Remark: there is one momentum integral for every interval (N in total), while there is one position integral for every intermediate position (N-1 in total). In the continuum limit, N\rightarrow\infty, this last expression approximates an integral over functions p(t),q(t), and we write it as

    \[ \displaystyle{\boxed{K\equiv \int \mathcal{D}p(t)\mathcal{D}q(t)e^{i\int _0^Tdt S(p,q)}}}\]

where

    \[ S(q,p)=p\dot{q}-H(p,q)\]

is the classical action in phase space! This phase integral is over all the classical functions p(t),q(t), with q(t) such as q(0)=q and q(T)=q'.

Now we are ready for two simple and frequently found examples of path integrals, very important indeed in practical applications!

Examples of path integrals

The first easy example is the free particle propagator and its path integral. The propagator K(q',q;T,0) will be related to the hamiltonian

H=\dfrac{p^2}{2m}

We can compute the propagator with simple QM:

K=\langle q'\vert e^{-iHT}\vert q\rangle

thus

\displaystyle{\langle q'\vert e^{-iT\frac{p^2}{2m}}\int\dfrac{dp}{2\pi}\vert p\rangle\langle p\vert q\rangle=\int\dfrac{dp}{2\pi}e^{-i\frac{Tp^2}{2m}}\langle q'\vert p\rangle\langle p\vert q\rangle}

so

\displaystyle{\boxed{K=\int\dfrac{dp}{2m}e^{-i\frac{Tp^2}{2m}+i(q'-q)p}}}

The integral is a gaussian, so we can integrate it inmediately, to obtain

\boxed{K=\left(\dfrac{m}{2\pi iT}\right)^{1/2}e^{im\frac{(q'-q)^2}{2T}}=\left(\dfrac{m}{2\pi iT}\right)^{1/2}e^{iS_c}}

and where the classical action S_c for a given time T reads

S_c=\dfrac{1}{2}mv^2T=\dfrac{p^2T}{2m}=\dfrac{1}{2}m\dfrac{(q'-q)^2}{T}

Of course, this approach is the normal QM procedure. We can get the same result using our brand new path integral approach!!! The configuration space PI reads:

\displaystyle{K=\lim_{N\rightarrow\infty}\left(\dfrac{m}{2\pi i\delta}\right)^{N/2}\int\prod_{j=1}^{N-1}dq_je^{im\frac{\delta}{2}\sum_{j=0}^{N-1}\left(\frac{q_{j+1}-q_{j}}{\delta}\right)^2}}

or expanding the last sum

\displaystyle{K=\lim_{N\rightarrow\infty}\left(\dfrac{m}{2\pi i\delta}\right)^{N/2}\int\prod_{j=1}^{N-1}dq_je^{im\frac{1}{2\delta}\left((q_N-q_{N-1})^2+\ldots+(q_1-q_0)^2\right)}}

where the initial and final points are q_0=q_i=q and q'=q_f=q_N. The integrals are all gaussian and can be evaluated exactly. The only issue is that they are “coupled” and it complicates the integration. You can do it as an interesting exercise. The final result is

\displaystyle{K=\lim_{N\rightarrow\infty}\left(\dfrac{m}{2\pi i\delta}\right)^{N/2}\dfrac{1}{\sqrt{N}}\left(\dfrac{2\pi i \delta}{m}\right)^{(N-1)/2}e^{im(q'-q)^2/2N\delta}=\lim_{N\rightarrow\infty}\left(\dfrac{m}{2\pi iN\delta}\right)^{1/2}e^{im(q'-q)^2/2\delta N}}

and taking into account that the total time is T=N\delta we write and check the previously obtained result

K=\left(\dfrac{m}{2\pi i T}\right)^{1/2}e^{iS_c}

with S_c=S(q_c) as the classical action as before. If we restore the Planck constant, the free particle propagator reads

\boxed{K=\left(\dfrac{m}{2\pi i\hbar T}\right)^{1/2}e^{iS_c/\hbar}}

and the propagator separates into two factors, one of which is the “classical phase” and the remaining one is a weird term depending on the square root of the mass and the time T.

The second example is the calculation of the path integral and propagator for the harmonic oscillator potential. We have:

\displaystyle{K(q',q;T,0)=\int \mathcal{D}q(t)e^{iS_c}}

The harmonic oscillator action reads

\displaystyle{S(q)=\int_0^Tdt\left(\dfrac{1}{2}m\dot{q}^2-\dfrac{1}{2}m\omega^2q^2\right)}

The paths are integrated out with q(0)=q=q_i and q(T)=q'=q_f. To perform the integral, suppose we do know the solution of the classical problem, q_c. That is:

\ddot{q_c}+\omega^2q_c=0

q_c(0)=q, q_c(T)=q'

so

q(t)=q_c(t)+y(t)

and we change variables in the path integral to

y(t)=q(t)-q_c(t)

Integrating out over all deviations from the classical path is equivalent to integrating out all the possible paths.

Remark: The jacobian of the above transformation is 1.

Remark(II): Since the classical solution q_c(t) obeys the correct boundary conditions, the paths y(t) over which we integrate go from

y(0)=0 to y(T)=0

The action can be written in terms of the path q_c+y=q as a power series in y(t):

\displaystyle{S=\int_0^Tdt\left(\dfrac{1}{2}m\dot{q_c}^2-\dfrac{1}{2}m\omega^2q_c^2\right)+\int_0^Tdt\left(\dfrac{1}{2}m\dot{y}^2-\dfrac{1}{2}m\omega^2y^2\right)+\mbox{linear terms in y}}

The terms linear in y vanish by construction. Then, we write

S(q_c+y)=S(q_c)+S(y)

By direct substitution, we get

\displaystyle{K(q',q;T,0)=e^{iS_c}\int \mathcal{D}y(t)e^{iS(y(t))}}

The path integral separates into two pieces:

1st. The term of the action with the classical action.

2nd. The path integral over the deviations from the classical path. This second factor is independent of the initial and final points!

We can calculate the path integral in position space, as we did above for the harmonic oscillator or, alternatively, we can compute it in the Fourier space, writing

\displaystyle{y(t)=\sum_ka_k\sin (k\pi t/T)}

and integrating over the coefficients a_k, the result will be

\boxed{K(q',q;T,0)=\left(\dfrac{m\omega}{2\pi i\sin\omega T}\right)^{1/2}e^{iS_c}}

Remark: For tiny frequencies, the harmonic oscillator propagator becomes the free particle propagator.

Remark(II): The classical action can be evaluated in a straightforward fashion, to get

S(q_c)=\dfrac{m\omega}{2\sin\omega T}\left((q'^2+q^2)\cos\omega T-2q'q\right)

A final comment follows: the path integral for ANY quadratic action can be evaluated EXACTLY, essentially performing a path integral with the aid of gaussian integrals!

Statistical Mechanics with Path Integrals

The central object in Statistical Mechanics is the partition function Z. Indeed, there is a beautiful correspondence between Statistical Mechanics and the path integral approach to QM by realizing the link between the partition function and the path integral (or the propagator). The classical partition function reads

\displaystyle{Z=\sum_je^{-\beta E_j}}

and where \beta=1/k_BT is the reciprocal temperature (in natural units) and E_j are the energies of the jth levels of the system, \vert j\rangle. We can write,

(1)   \begin{equation*} \displaystyle{Z=\sum_j\langle j\vert e^{-\beta H}\vert j\rangle=\mbox{Tr}e^{\beta H}}\end{equation*}

But, the propagator definition was

K(q',q;T,0)=\langle q'\vert e^{-iHT}\vert q\rangle

Suppose we consider T to be a COMPLEX parameter (it sounds a bit weird since time is real, not complex) and consider it to be a pure imaginary number. In particular, we write a “Wick” rotation of time as the temperature!

\boxed{T=-i\beta=\dfrac{\beta}{i}=\dfrac{1}{ik_B t}}

where t is the absolute temperature (be aware, here t is NOT time but temperature and, though, they are closely related) and it is a REAL number. Therefore,

\displaystyle{K(q',q;-i\beta,0)=\langle q'\vert e^{-iH(-i\beta)}\vert q\rangle=\langle q'\vert e^{-\beta H}\sum_j\vert j\rangle\langle j\vert q\rangle=\sum_je^{-\beta E_j}\langle q'\vert j\rangle\langle j\vert q\rangle=\sum_je^{-\beta E_j}\langle j\vert q\rangle\langle q'\vert j\rangle}

Inserting q'=q and integrating over the position variable q, we obtain

\displaystyle{\int dq K(q,q;-i\beta,0)=\sum_j e^{-\beta E_j}\langle j\vert \int dq\vert q\rangle\langle q\vert j\rangle=Z}

And thus, we have the formal identity and link between the path integral (the propagator) and the partition function of statistical mechanics:

\displaystyle{\boxed{Z=\int dq K(q,q;-i\beta,0)}}

Put it into (equivalent) words:

1st. The propagator evaluated at negative imaginary time (a Wick rotated time) IS the partition function. 

2nd. Temperature is a Wick rotation of time or time is a Wick rotated temperature.

3rd. The propagator in diagonal position representation IS the partition function.

Fascinating! We can easily work out an elementary example with the harmonic oscillator above. Recall the PI for it:

K(q',q;T,0)=\left(\dfrac{m\omega}{2\pi i\sin\omega T}\right)^{1/2}\exp\left(i\dfrac{m\omega}{2\sin\omega T}\left((q'-q)^2\cos\omega T-2q'q\right)\right)

Write the Wick rotation T=-i\beta and q'=q (diagonal representation of the partition function/propagator):

K(q,q;-i\beta,0)=\left(\dfrac{m\omega}{2\pi\sinh(\beta\omega)}\right)^{1/2}\exp\left(-\dfrac{m\omega q^2}{\sinh(\beta\omega)}(\cosh(\beta\omega)-1)\right)

This propagator gives a partition function

\displaystyle{Z=\int dqK(q,q;-i\beta,0)=\left(\dfrac{m\omega}{2\pi\sinh(\beta\omega)}\right)^{1/2}\sqrt{\dfrac{\pi}{\frac{m\omega}{\sinh(\beta\omega)}(\cosh(\beta\omega)-1)}}=\sqrt{\dfrac{1}{2(\cosh(\beta\omega)-1)}}}

and then

\displaystyle{Z=\dfrac{1}{e^{\beta\omega/2}\left(1-e^{-\beta\omega}\right)}}

But this is precisely the classical (quantum) result

\displaystyle{Z=\dfrac{e^{-\beta\omega/2}}{1-e^{-\beta\omega}}=\sum_{j=0}^{\infty}e^{-\beta (j+\frac{1}{2})\omega}=\sum_{j=0}^{\infty} e^{-\beta\hbar\omega (j+\frac{1}{2})}}

and where we reinserted the Planck constant in the last step! Path integrals are disguised partition functions and vice versa!

We can rewrite the partition function in terms of a path integral. In ordinary (real) time,

\displaystyle{K(q',q;T,0)=\int\mathcal{D}q(t)\exp\left(i\int_0^Tdt\left(\dfrac{m\dot{q}^2}{2}-V(q)\right)\right)}

and where the integral is taken over all paths from (q,0) to (q',T). In the diagonal representation, q'=q and T=-i\beta is the Wick rotated time/temperature, so

\displaystyle{K(q,q;-i\beta,0)=\int\mathcal{D}q\exp\left(i\int_0^{-i\beta}dt\left(\dfrac{m\dot{q}^2}{2}-V(q)\right)\right)}

and the integral in now along the negative imaginary time/temperature axis! Let us define a real variable for this integration, say \tau=it. Here, \tau is called the imaginary time, since when the time t is IMAGINARY, then \tau is REAL. Then, the integral over the IMAGINARY time \tau is along the real axis t: when t runs from 0 to -i\beta, then \tau runs from 0 to \beta. We can write q as a function of the variable \tau: q(t)\rightarrow q(\tau)=q(it). It defines the so-called imaginary time formalism. Then,

    \[ \dot{q}=i\dfrac{dq}{d\tau}=p\]

Note the similarity of the above last equation to the time dependent Schrödinger equation and the classical Hamiltonian Dynamics! Then, the propagator becomes

    \[ \displaystyle{K(q,q;-i\beta,0)=\int\mathcal{D}q\exp\left(-\int_0^\beta d\tau\left(\dfrac{m}{2}\left(\dfrac{dq}{d\tau}\right)^2+V(q)\right)\right)}\]

The integral is taken over all function q(\tau) such that

q(\beta)=q(0)=q

This final result is a Wick rotated “imaginary time” or “euclidean” path integral, defined by associating to any path an amplitude (statistical weight) e^{-S_E}, where S_E is the so-called euclidean action, obtained from the usual (or Minkovskian, spacetime) action by changing the sign of the potential energy term.

In summary, we have got two main formulae:

1st. The Minkovski or real time propagator, in which we have an oscillating integral of the classical action

\displaystyle{K(q',q;T,0)=\int\mathcal{D}q(t)\exp\left(i\int_0^Tdt\left(m\dfrac{\dot{q}^2}{2}-V(q)\right)\right)}

2nd. The Euclidean or imaginary time propagator, in which we have an exponential integral of the classical action

\displaystyle{K(q,q;-i\beta,0)=\int\mathcal{D}q(\tau)\exp\left(-\int_0^\beta d\tau \left(\dfrac{m}{2}\left(\dfrac{dq}{d\tau}\right)^2+V(q)\right)\right)}

The integral here is over all functions q(\tau) such as q(0)=q(\tau)=q

We can use one of another formula, but usually euclidean time are preferred in most of the applications due to convergence features. We could perform all the calculations in a Euclidean theory and analytically continue the results to the real time when computing the physical quantities as well!

From the partition function, the statistical partition function, it is possible to extract the ground state energy. In fact, from the definition of Z

    \[ \displaystyle{Z(\beta)=\sum_j e^{-\beta E_j}=e^{-\beta E_0}(1+e^{-\beta \Delta E_j})}\]

and

    \[ \Delta E_j\equiv E_j-E_0\]

We can see that the contribution of every state decreases exponentially with \beta, the temperature (reciprocal of). However, the contribution of the energy of the ground state decreases less slowly than any other state contribution. Thus, in the limit of large \beta, i.e. low temperature, the ground state contribution will dominate. We arrive at

\displaystyle{E_0=-\lim_{\beta\rightarrow\infty}=\dfrac{1}{\beta}\log Z}

We can do the same trickery in an alternative way. Let us look at the euclidean time propagator from q=0 to q'=0:

K(0,0;\beta,0)=K_E(0,0;\beta,0)=\langle q'=0\vert e^{-\beta H}\vert q=0\rangle=\langle 0\vert e^{-\beta H}\vert 0\rangle

We can insert a complete set of eigenstates of H:

\displaystyle{K_E(0,0;\beta,0)=\langle 0\vert e^{-\beta H}\sum_j\vert j\rangle\langle j\vert 0\rangle=\sum_j e^{-\beta E_j}\vert \psi_j(0)\vert^2}

or

\displaystyle{K_E(0,0;\beta,0)=e^{-\beta E_0}\vert \psi_0(0)\vert^2\left(1+\sum_{j=1}^\infty e^{-\beta\Delta E_j}\vert\psi_j(0)/\psi_0(0)\vert^2\right)}

and where \psi_j(0) are the wave functions of H. Again, the ground state dominates at low temperatures (\beta\rightarrow\infty), and we get

\displaystyle{E_0=-\lim_{\beta\rightarrow\infty}\dfrac{1}{\beta}\log (K_E(0,0;\beta,0)-\vert\psi_0(0)\vert^2)=-\lim_{\beta\rightarrow\infty}\dfrac{1}{\beta}\log K_E(0,0;\beta,0)}

Namely, in the low temperature regime the difference between the terms \beta^{-1}\log Z and \beta^{-1}\log K_E goes to zero. And this result is important…

Density matrix

The density matrix is an important object. It has a well known expression in quantum statistical mechanics:

\boxed{\rho=\dfrac{e^{-\beta H}}{\mbox{Tr}(e^{-\beta H})}}

Using the previous results, we have

\rho (q',q)=\dfrac{1}{Z}K(q',q;-i\beta,0)

Indeed, every observable can be related to the partition function, the path integral, and the density matrix with these definitions! The expectation value for some operator \hat{O} is given by

\boxed{\langle \hat{O}\rangle=\mbox{Tr}(\hat{O}\rho)}

Using the path integral, we can write

\displaystyle{\boxed{\langle \hat{O}\rangle=\dfrac{\int \mathcal{D}q \hat{O}e^{-S_E}}{\int\mathcal{D}q e^{-S_E}}}}

Correlators in QM and QFT

The position operator in the usual Heisenberg representation q(t) is defined in terms of the Schrödinger operator picture by means of the identity

q_H(t)=e^{iHt}q_Se^{-iHt}

The time eigenstates read

q(t)\vert q,t\rangle=q\vert q,t\rangle

and its relation with the time-independent eigenstates is

\vert q,t\rangle=e^{iHt}\vert q\rangle

We can express the path integral with these states:

\displaystyle{K=\langle q'\vert e^{-iHT}\vert q\rangle=\langle q',T\vert q,0\rangle=\int\mathcal{D}qe^{iS}}

We can define the time ordering of operators as follows. Let us have two time dependent operators A(t), B(t). The time ordered operator (TOO) is:

TA(t_1)B(t_2)=\begin{cases}A(t_1)B(t_2),t_1>t_2\\ B(t_2)A(t_1),t_1<t_2\end{cases}

We also define the n-order Green function, or n-point correlator function as the following expectation value in QM:

\boxed{G^{(n)}(t_1,t_2,\ldots,t_n)=\langle 0,t_f=\infty \vert Tq(t_1)q(t_2)\cdots q(t_n)\vert 0,t_i=-\infty\rangle}

An important case is the so-called 2-point (correlation) function. The two point function G(t_1,t_2) has a very special importance in QM. Let us now calculate the 2-point function G(t_1,t_2) via the path integral formalism:

G(t_1,t_2)=\langle 0,\infty\vert Tq(t_1)q(t_2)\vert 0,-\infty\rangle

Firstly, we will calculate

\langle q',T\vert Tq(t_1)q(t_2)\vert q,-T\rangle

and then, we will take the limits q'=q=0 and T\rightarrow\infty. Suppose that t_1>t_2, then

\displaystyle{\langle q',T\vert Tq(t_1)q(t_2)\vert q,-T\rangle=\int dq_1dq_2\langle q',T\vert q_1,t_1\rangle\langle q_1,t_1\vert q_2,t_2\rangle\langle q_2,t_2\vert q,0\rangle}

Every matrix element is itself a path integral:

\displaystyle{\langle q',T\vert Tq(t_1)q(t_2)\vert q,-T\rangle=\int dq_1 dq_2\int_{q_1,t_1}^{q',T}\mathcal{D}q e^{iS}\int_{q_2,t_2}^{q_1,t_1}\mathcal{D}qe^{iS}\int_{q,0}^{q_2,t_2}\mathcal{D}qe^{iS}}

This last expression have 3 different pieces:

1st. A path integral form the initial position q to an arbitrary position q_2.

2nd. A path integral from q_2 to a second arbitrary position q_1.

3rd. A path integral from an arbitrary position q_1 to the final position at q'.

We are integrating over all paths from q to q’, subject to the restriction that the paths pass through the intermediate points at q_1 and q_2. We then integrate over the two arbitrary positions, so that in fact we are integrating over ALL paths: we can combine these three path integrals plus the integrations over q_1 and q_2 into one single path integral. The factos q_1,q_2 in the above integral can be incorporated into this simple integral by including a factor q_1,q_2=q_1q_2 into the path integral. Therefore, we read

\displaystyle{\langle q',T\vert q(t_1)q(t_2)\vert q,-T\rangle=\int_{q,-T}^{q',+T}\mathcal{D}qq(t_1)q(t_2)e^{iS},t_1>t_2}

Mimicking this calculation, in the case of t_1<t_2:

\displaystyle{\langle q',T\vert q(t_2)q(t_1)\vert q-T\rangle=\int_{q,-T}^{q',+T}\mathcal{D}qq(t_1)q(t_2)e^{iS}}

That is, the same expression is valid in both cases! The path integral does the time ordering automatically!!!! Taking now the limit T\rightarrow\infty, we get

\displaystyle{\boxed{G^{(2)}(t_1,t_2)=\langle 0\vert Tq(t_1)q(t_2)\vert -0\rangle=\int_{-0,0}\mathcal{D}qq(t_1)q(t_2)e^{iS}}}

For the n-point correlation functions, we obtain similar results:

\displaystyle{\boxed{G^{(n)}(t_1,t_2,\ldots,t_n)=\langle +0\vert Tq(t_1)q(t_2)\cdots q(t_n)\vert -0\rangle=\int_{-0,+0}\mathcal{D}qq(t_1)q(t_2)\cdots q(t_n)e^{iS}}}

 The generating functional Z(J)

In order to define a generating function Z(J), we have to build up some preliminary definitions. The general amplitude

A=\langle x\vert U(t,t_0)\vert x_0\rangle\equiv \langle x,t\vert x_0,t_0\rangle

The vacuum amplitude or 0-0 vacuum persistence amplitude is:

    \[ \boxed{\langle +0\vert -0\rangle=\langle 0,t=+\infty\vert 0,t=-\infty\rangle}\]

with

    \[ \displaystyle{\langle +0\vert -0\rangle=\int\mathcal{D}qe^{iS}}\]

where the integration is calculated over all the closed paths such that

q(t_0=-\infty)=q(t=+\infty)=0

We will see how this amplitude is related to the ground state of survival (or persistence) amplitude. Firstly, we add a new “source” to the action of an arbitrary system. Thus, we modify the action with an extra term/piece:

\displaystyle{S(q,J)=S(q)+\int dtJ(t)q(t)}

and where J(t) is an arbitrary function such as

J(t\rightarrow\pm\infty)=0

The generating functional Z(J) is defined now to be the quantity

\displaystyle{\boxed{Z(J)=\dfrac{\int\mathcal{D}qe^{iS+i\int dt J(t)q(t)}}{\int\mathcal{D}qe^{iS}}=\dfrac{\langle +0\vert -0\rangle_J}{\langle +0\vert -0\rangle_{J=0}}}}

with q(t_0=-\infty)=q(\infty)=0 as before. The functional derivatives of the generating functional Z(J) are very interesting objects:

\dfrac{1}{i}\dfrac{\delta Z(J)}{\delta J(t_1)}

evaluated at J=0 provides the one-point correlation function! The proof is simple:

\left(\dfrac{1}{i}\dfrac{\delta Z(J)}{\delta J(t_1)}\right)\vert_{J=0}=\left(\dfrac{\int \mathcal{D}q(t_1)e^{iS+i\int dtJ(t)q(t)}}{\int\mathcal{D}qe^{iS}}\right)\vert_{J=0}

or equivalently

\left(-i\dfrac{\delta Z(J)}{\delta J(t_1)}\right)\vert_{J=0}=\dfrac{\int\mathcal{D}qq(t_1)e^{iS}}{\int\mathcal{D}qe^{iS}}=\langle +0\vert q(t_1)\vert -0\rangle=G^{(1)}(t_1)

Iterating this process, we obtain a path integral with several q’s in the numerator, i.e., we get the n-point (correlation) function!

\boxed{\left((-i)^n\dfrac{\delta}{\delta J(t_1)}\cdots\dfrac{\delta}{\delta J(t_{n-1})}\dfrac{\delta Z(J)}{\delta J(t_n)}\right)\vert_{J=0}=\dfrac{\int\mathcal{D}qq(t_1)q(t_2)\cdots q(t_n)e^{iS}}{\int\mathcal{D}qe^{iS}}=\langle +0\vert Tq(t_1)q(t_2)\cdots q(t_n)\vert -0\rangle=G^{(n)}(t_1,t_2,\ldots,t_n)}

Indeed, using the functional derivatives and expanding the generating functional Z(J) in Taylor-like series with these derivatives in a series on J, we get the formal series:

Z(J)=Z(0)+Z(1)+Z(2)+\cdots

where

Z(0)=\langle +0\vert -0\rangle

Z(1)=\int dyG^{(1)}(y)J(y)

Z(2)=\dfrac{1}{2!}\int dy_0dy_1G^{(2)}(y_0,y_1)J(y_0)J(y_1)

and so on! You can understand now why Z(J) is called the generating functional: its expansion provides the n-point correlation functions!

We can compute the generating functional of the harmonic oscillator as a working example. Suppose that initially the action of this system reads

\displaystyle{\int dt\left(\dfrac{1}{2}m\dot{q}^2-\dfrac{1}{2}m\omega^2q^2\right)}

Then, the corresponding denominator, N(0), is the propagator K we computed for the harmonic oscillator before! For the numerator in the generating functional, we write

\displaystyle{N(J)=\mathcal{D}e^{iS+i\int dt J(t)q(t)}}

By definition, we also get

\displaystyle{N(J)=\int\mathcal{D}q\exp\left(i\int dt\left(\dfrac{1}{2}m\dot{q}^2-\dfrac{1}{2}m\omega^2q^2+Jq\right)\right)}

Now, we use the same mathematical magic we did for the unforced harmonic oscillator, defining

q(t)=q_c(t)+y(t)

and we integrate over all the possible paths and the classical solution q_c satisfies

m\ddot{q_c}+m\omega^2q_c=J(t)

q(t_0=-\infty)=q(t=\infty)=0

and then

\displaystyle{\int dt\left(\dfrac{1}{2}m\dot{q}^2-\dfrac{1}{2}m\omega^2q^2+J(t)q(t)\right)=\int dt\left(\dfrac{1}{2}m\dot{q_c}^2-\dfrac{1}{2}m\omega^2q_c^2+J(t)q_c(t)\right)+\int dt\left(\dfrac{1}{2}m\dot{y}^2+\dfrac{1}{2}m\omega^2y^2\right)+\mbox{linear terms in y}}

The linear terms in y vanish since q_c satisfies the equation of motion. Therefore, the path integral reduces to

\displaystyle{N(J)=e^{iS_c}\int\mathcal{D}y\exp\left(i\int dt\left(\dfrac{1}{2}m\dot{y}^2-\dfrac{1}{2}m\omega^2y^2\right)\right)}

And now, the path integral over y is a constant, independent from J, that we will call C. We write

N(J)=Ce^{iS_c(q_c,J)}

where

\displaystyle{S(q_c,J)=\int dt \left(\dfrac{1}{2}m\dot{q_c}^2-\dfrac{1}{2}m\omega^2q_c^2+Jq_c\right)}

We can compute this last integral with the aid of the equation of motion and integrating the first term by parts:

\displaystyle{\int dt (-1)q_c(t)\dfrac{1}{2}\left(m\dfrac{d^2}{dt^2}+m\omega^2\right)q_c(t)+J(t)q_c(t)=\int dt (-1)\dfrac{1}{2}q_c(t)J(t)+J(t)q_c(t)}

and it yields

\displaystyle{S(q_c,J)=\dfrac{1}{2}\int dt J(t)q_c(t)}

We can also determine the classical path and define the Green function with this new formalism. The classical equation of motion is

\ddot{q_c}+\omega^2q_c=J(t)

The classical path can be written in terms of the Green functions defined by the formal identities

\left(\dfrac{d^2}{dt^2}+\omega^2\right)iG(t',t)=\delta (t'-t)

G(t',t=-\infty)=G(t',t=+\infty)=0

and

q_c(t)=-i\int dt'G(t',t)J(t')

We can solve these equations for the Green function by going into the momentum space representation by Fourier transformations. The result is

\displaystyle{G(t,t')=G(t-t')=\int\dfrac{dk}{2\pi}\dfrac{i}{k^2-\omega^2}e^{-ik(t-t')}}

However, there are poles on the axis of integration. The Green function is ambiguous unless we give it a pole prescription, i.e., a boundary condition (namely, we have to avoid the pole singularity). We require that G go to zero as T\rightarrow\infty. The correct pole prescription then turns out to be

\displaystyle{G(t-t')=\int\dfrac{dk}{2\pi}\dfrac{i}{k^2-\omega^2+i\varepsilon}e^{-ik(t-t')}}

Note that in the case of equal times, the integral is FINITE and the Green function simplifies to

G(t,t)=G(0)=\dfrac{1}{2\omega}

We can now write

\displaystyle{N(J)=C\exp\left(\dfrac{1}{2}\int dt'dtJ(t)G(t,t')J(t')\right)}

Dividing by the denominator merely cancels the constant C factor, and thus

Z(J)=\dfrac{N(J)}{N(0)}

and the final result will be

\displaystyle{Z(J)=\exp\left(\dfrac{1}{2}dtdt'J(t)G(t,t')J(t')\right)}

and where G is the Green function defined above!

Now we can study some general features or properties. Let me take any functional K(J):

\displaystyle{K(J)=\int\mathcal{D}qe^{iS(q)}=\exp\left(i\left(S(q)+\int dtJ(t)q(t)\right)\right)}

We can compute the first functional derivative for any time t_1

\displaystyle{\dfrac{\delta K(J)}{\delta J(t_1)}=\int\mathcal{D}qq(t_1)\exp\left(i\left(S(q)+\int dtJ(t)q(t)\right)\right)}

The second derivative can be also calculated further

\displaystyle{\dfrac{\delta^2K(J)}{\delta J(t_1)\delta J(t_2)}=\int\mathcal{D}qq(t_1)q(t_2)e^{iS(q,t)}}

Indeed, we can generalize this result to an arbitrary functional F in the following way

\displaystyle{F\left(\dfrac{\delta K(J)}{\delta J}\right)=\int\mathcal{D}qF(q)e^{iS(q,J)}}

This expression can be proved correct in two steps:

1st. Put F(\delta/\delta J) inside the PI, since it does not depend on the path.

2nd. Expand F in a functional Taylor series, and every functional derivative with respect to J provides a q term in front of the exponential term.

Remark: As an example, take

\displaystyle{F(q)=\exp\left(i\left(\int dtV(q)\right)\right)}

We get in this case

\displaystyle{e^{i\int dt V(\delta/\delta J(t))}K(J)=\int\mathcal{D}qe^{i\int dtV(q)}e^{iS(q,J)}=\int\mathcal{D}qe^{iS(q,J)+i\int dtV(q)}}

The previous expression is the starting point of perturbation theory in QM! Of course, we are making perturbation theory with the path integral formalism, not in the usual way but it looks pretty similar and it is done in a completely covariant way. That is the hidden power of path integrals. Let us expand the exponential in the right-handed side above and write the expression as follows

\displaystyle{e^{i\lambda\int dtV\left(\delta/\delta J(t)\right)}=1+i\lambda\int dtV\left(\dfrac{\delta}{\delta J(t)}\right)+\dfrac{1}{2}\left(i\lambda\int dtV\left(\dfrac{\delta}{\delta J(t)}\right)\right)^2+\cdots}

or

\displaystyle{e^{i\lambda\int dtV\left(\delta/\delta J(t)\right)}=1+i\lambda\int dtV\left(\dfrac{\delta}{\delta J(t)}\right)+\dfrac{1}{2}(i\lambda)^2\int dt_1V\left(\dfrac{\delta}{\delta J(t_1)}\right)\int dt_2V\left(\dfrac{\delta}{\delta J(t_2)}\right)+\cdots}

or

\displaystyle{e^{i\lambda\int dtV\left(\delta/\delta J(t)\right)}=1+i\lambda\int dtV\left(\dfrac{\delta}{\delta J(t)}\right)+\dfrac{1}{2}(i\lambda)^2\int dt_1dt_2V\left(\dfrac{\delta}{\delta J(t_1)}\right)V\left(\dfrac{\delta}{\delta J(t_2)}\right)+\cdots}

This is the analogue form of the so-called Dyson series in the PI formalism! For small values of the parameters \lambda, we get an approximation to the PI truncating the series at any desired order. At first order in \lambda we write

\displaystyle{\int\mathcal{D}qe^{iS(q,J)+i\lambda\int dtV(q)}\simeq\left(1+i\lambda\int dtV\left(\dfrac{\delta}{\delta J(t)}\right)\right)K(J)=K(J)+i\lambda\int dt V\left(\dfrac{\delta}{\delta J(t)}\right)K(J)}

If we plug the condition J=0 at the end of our calculation, it yields

\displaystyle{\int\mathcal{D}qe^{iS(q,J)+i\lambda\int dtV(q)}\simeq\int\mathcal{D}qe^{iS(q)}+i\lambda\int dtV\left(\dfrac{\delta}{\delta J(t)}\right)K(J)\vert_{J=0}}

A very important popular (VIP) example is provided by the quartic interaction potential. It provides a functional term

\displaystyle{F(q)=\exp\left(i\int d\tau\dfrac{\lambda}{4!}q^4\right)}

The \lambda \phi^4 potential term is very important in the Standard Model and Quantum Mechanics. We obtain

\displaystyle{\int\mathcal{D}qe^{iS(q)+i\int dt\lambda\frac{1}{4!}q^4}=\int\mathcal{D}qe^{iS(q)}+i\dfrac{\lambda}{4!}\int dt\left(\dfrac{\delta}{\delta J(t)}\right)^4K(J)\vert_{J=0}+o(\lambda^2)}

and where

\displaystyle{\boxed{K(J)=\int\mathcal{D}qe^{iS(q)+i\int dtJ(t)q(t)}}}

Have you enjoyed this post? See you in my last introductory path integral TSOR post…

P.S.: A special post, the 150th, is coming soon! Stay tuned!

View ratings
Rate this article

Comments

LOG#148. Path integral (III). — 2 Comments

  1. First, nice blog!

    Second, I’ve watched a conference held by Prof. Unruh a few months ago, where he claimed the experimental proof of his effect in the vibration field on the surface of a flowing liquid, i.e. an “acoustic black hole”. He argued that his and the Hawking effect are in fact one and the same, and that there’s no real “quantum” origin behind the black hole thermal radiation. He also dismissed the whole “virtual particle pair creation at the event horizon” as something that he considers unphysical nonsense. I don’t remember the details though as GR and quantum gravity are not my field, I guess his papers probably explain it better. What do you think about this? In a sense, considering this to be a natural effect that takes place in all fields without any intervention from quantum effects, and that it arises in acoustic fields as well as in extreme conditions like black holes would be a reassuringly elegant solution to the issues at hand. In another, it seems weird that something so trivial as waves on a liquid surface can arise from the same reasons that cause something as crazy as black holes to irradiate (and potentially evaporate), when we know that at the event horizon there might as well be effects we don’t understand at play, given how we can’t reconcile GR and quantum mechanics.
    BTW, the paper by Unruh describing the experiment can be found here:

    http://arxiv.org/abs/1008.1911

    Reply ↓

    • Well, as I have told you in facebook. I am skeptical about what Unruh is defending now. Who knows? Who nose? Unless experimental evidence of Unruh effect and other quantum gravitational experiments involving the vacuum appears, we have only theories and conjectures…:) Best wishes!

      Reply ↓

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.