# Path Integrals: From Quantum Mechanics to Photonics

Charles W. Robson  Photonics Laboratory, Tampere University, P.O. Box 692,
FI-33014, Tampere, Finland  Yaraslau Tamashevich  Photonics Laboratory,
Tampere University, P.O. Box 692, FI-33014, Tampere, Finland  Tapio T. Rantala
Department of Physics, Tampere University, P.O. Box 692, FI-33014, Tampere,
Finland  Marco Ornigotti<EMAIL_ADDRESS>Photonics Laboratory,
Tampere University, P.O. Box 692, FI-33014, Tampere, Finland

###### Abstract

The path integral formulation of quantum mechanics, __i.e.__ , the idea that
the evolution of a quantum system is determined as a sum over all the possible
trajectories that would take the system from the initial to its final state of
its dynamical evolution, is perhaps the most elegant and universal framework
developed in theoretical physics, second only to the Standard Model of
particle physics. In this tutorial, we retrace the steps that led to the
creation of such a remarkable framework, discuss its foundations, and present
some of the classical examples of problems that can be solved using the path
integral formalism, as a way to introduce the readers to the topic, and help
them get familiar with the formalism. Then, we focus our attention on the use
of path integrals in optics and photonics, and discuss in detail how they have
been used in the past to approach several problems, ranging from the
propagation of light in inhomogeneous media, to parametric amplification, and
quantum nonlinear optics in arbitrary media. To complement this, we also
briefly present the Path Integral Monte Carlo (PIMC) method, as a valuable
computational resource for condensed matter physics, and discuss its potential
applications and advantages if used in photonics.

††preprint: AIP/123-QED

##  I Introduction

The path integral formulation of quantum mechanics, developed in the mid-
twentieth century, is not only a remarkable synthesis of several of the core
ideas of theoretical physics, but it is also a powerful computational
technique for the analysis of a huge variety of physical systems in very
different contexts, such as quantum mechanics Feynman (1948); Feynman and
Hibbs (1965); Zinn-Justin (2005); (2005) (ed.); Kleinert (2007); Schulman
(2005), quantum field theory Fai (2019); Das (2006); Rivers (1988); Maggiore
(2005), gauge field theory Fadeev and Slavnov (1991); Frampton (2008), black
hole physics Schulman (2005), quantum gravity Hawking (1978); Birrell and
Davies (1982), string theory Polchinski (2005), topology Schulman (2005);
Mourao (2004)condensed matter physics Pollock and Ceperley (1984); Ceperley,
Simmons, and Blasdell (1996); Ceperley (1995, 1996); Pierleoni  _et al._
(1994); Militzer and Ceperley (2001); Gorelov  _et al._ (2020), statistical
mechanics Feynman (1998); Bhagwat, Khandekar, and Lawande (1993); Feynman and
Hibbs (1965); Schulman (2005), polymer physics Kleinert (2007); Schulman
(2005), financial markets Kleinert (2007), optical communications Reznichenko,
Terekhov, and Turitsyn (2017), atomic physics Zaheer and Zubairy (1988);
Chernyak and Mukamel (1994), spectroscopy Tanimura and Mukamel (1993), light
propagation in turbid media Perelman  _et al._ (1994) and classical Gómez-
Reino and Liñares (1987); Braidotti  _et al._ (2018); Dimant and Levit (2010);
Gersten and Nitzan (1987); Wright and Garrison (1987); Vatsya (2005) and
quantum Hillery and Zubairy (1984, 1982); Difallah, Szameit, and Ornigotti
(2019) optics, amongst others.

The underlying concepts of the path integral approach are sometimes considered
difficult to grasp. Indeed, from a philosophical standpoint, its underpinnings
are extraordinary – in describing a system, all possible paths between its
initial and final state must be taken into account mathematically. This idea
of summing over all paths has previously been characterised ontologically as
_everything that can happen does happen_ Cox and Forshaw (2012).

As befits the subject, the history of path integrals is rather circuitous and
does not follow a straight line from Feynman’s seminal work on the subject,
published in 1948 Feynman (1948), to the present day. Already in the 1920s,
the mathematician Norbert Wiener had developed a method to treat Brownian
motion and diffusion using a technique of integrating over paths, with formal
similarities to Feynman’s eventual construction but in a purely classical
context Wiener (1924); Brush (1961); Chaichian and Demichev (2001). During the
early part of the following decade, the general ideas of de Broglie and
Schrödinger, that waves can be associated with particle dynamics Bacciagaluppi
and Valentini (2009), motivated Dirac to publish a paper in the _Physikalische
Zeitschrift der Sowjetunion_ (Physical Journal of the Soviet Union) setting
the stage for future developments by proposing the Lagrangian as a more
natural basis for a theory of quantum mechanics, as opposed to a Hamiltonian-
based method which he argued was less fundamental due to its non-relativistic
form Dirac (1933). This led him to propose the key idea that a quantum
mechanical transition amplitude for a particle is given by a phase factor
controlled by the action along that particle’s path Dirac (1947). This also
led him to state, that the classical path of a quantum system can be
interpreted as resulting from the constructive interference of all such paths.
In other terms, the action $S$ of a given system counts, _de-facto_ , the
number of waves of the path in units of Planck’s constant $h$, and, therefore,
for each path, the phase of the related wave 111It is not difficult, using an
electromagnetic analogy, to see that in this context, Planck’s constant $h$
plays the role of wavelength of the path. A plane wave propagating, say, along
the $x$-direction, in fact, is associated with a phase factor of the form
$\exp{(i\,2\pi\,x/\lambda)}$. The phase factor carried by each path i,
instead, goven by $\exp{(i\,2\pi\,S/h)}$. A simple comparison allows the
reader to identify $\lambda\rightarrow h$. is given by
$\exp{(i\,2\pi\,S/h)}=\exp{(iS/\hbar)}$. This allows the evaluation of the
interference pattern of the particle dynamics, and essentially represents what
is nowadays commonly understood as path integral Feynman and Hibbs (1965);
Kleinert (2007).

The pivotal step in the development of the theory occurred in the 1940s, when
Feynman formulated his version of quantum mechanics, a ‘third way’ following
the earlier, well-known Schrödinger and Heisenberg alternatives (2005) (ed.).
Based on integration over all paths between initial and final physical states,
with each path contributing an action-dependent phase as Dirac had proposed,
Feynman’s formulation culminated in his important 1948 paper entitled _Space-
time approach to non-relativistic quantum mechanics_ Feynman (1948).

From a conceptual point of view, one of the most interesting consequences of
path integrals is that they can provide a deep understanding of the relation
between quantum and classical mechanics, as the limit $h\rightarrow 0$ emerges
naturally from the formalism as the classical limit of the theory Zinn-Justin
(2005).

Although the concepts underlying the path integral formalism may at first
appear quite alien, they are deeply profound. Even if its utility were
limited, the beauty of the idea would still merit a wide audience. The fact
is, however, that the path integral has immense value as a practical
computational technique in the physical sciences, and it can be applied to
solve problems in many diverse areas, even beyond physics. Two examples
showing the remarkable breadth of its applicability are its use in quantum
gravity, where the “sum over histories” is interpreted as a sum over all
different spacetime configurations interpolating between an initial and final
state of the universe Hawking (1978), and in financial market modelling, where
the formalism has proven useful, since the time dependence of asset prices can
be represented by fluctuating paths Kleinert (2007).

The universality of the technique has allowed scientists to tackle many
problems, and gain tremendous physical insight into them. In statistical
physics, for example, path integrals conveyed the basic framework for the
first formulation of the renormalisation group transformation, and they are
largely employed to study systems with random distribution of impurities
Schulman (2005). In particle physics, they allow one to understand and
properly account for the presence of instantons Weinberg (2005). In quantum
field theory they provide the natural framework to quantise gauge fields
Fadeev and Slavnov (1991); Frampton (2008); Ammon and Erdmenger (2015). In
chemical, atomic, and nuclear physics, on the other hand, they have been
applied to various semiclassical schemes for scattering theory Schulman
(2005). Through path integrals, topological and geometrical features of
classical and quantum fields can be readily investigated and be used to create
novel forms of perturbative and nonperturbative analysis of fundamental
processes of Nature Baez and Muniain (1994); Polchinski (2005); Green,
Schwarz, and Witten (2012). In addition to that, path integrals allow one to
re-interpret established results, such as, for example, the BCS theory of
superconductivity Fai (2019); Tempere and Devreese (2012), from a novel, more
insightful, perspective.

Explicit, analytical solutions to problems formulated in terms of path
integrals, however, are scarce, and only available for very simple systems,
such as a free particle, or the ubiquitous harmonic oscillator Feynman and
Hibbs (1965). The complexity of the path integral formalism, in fact,
increases very rapidly to overwhelming levels of difficulties for many simple
problems. As an example of that, the simplest quantum system, _i.e._ , a
single hydrogen atom, required nearly 40 years to be fully solved in terms of
path integrals Duru and Kleinert (1982).

On the other hand, it is amidst complex and computationally challenging
problems that path integrals show their true potential, providing a simple,
insightful and intuitive perspective on the physical principles regulating
such processes. To do that, numerical techniques, such as Monte Carlo methods
Creutz and Freedman (1981); Morningstar (2007); Creutz (1983); Gattringer and
Lang (2010); Westbroek  _et al._ (2018); Ceperley (1995); Gorelov  _et al._
(2020), and the computational power of modern supercomputers are crucial to
their successful implementation.

Most of the practical applications to numerical simulations of quantum
particles with path integral approaches are systems with finite temperature in
thermodynamical equilibrium with one of the statistical ensembles Feynman
(1998). Finite temperature equilibrium involves dynamics of constituent
particles and exchange of energy with environment or heat bath. These are the
central factors in condensed matter physics with phase transitions,
conductivities and other processes related to interactions between constituent
particles. Typically, one assumes canonical ensembles, where both the number
of particles and the volume they occupy are kept constant at a given
temperature $T$, but other ensembles can be chosen where needed.

For a many-particle system in finite temperature, there is no wave function,
but the mixed state can be described with a density matrix, and it turns out
that it can be written in terms of path integrals in imaginary time Feynman
and Hibbs (1965); Feynman (1998). The expectation values of observables are
then evaluated by using the trace of the density matrix, which in space basis
means finite closed loop paths. Thus, in terms of path integrals each of the
particles propagate from a position in space in imaginary time back to the
same position, the time period being inversely proportional to the
temperature. Then, with Metropolis Monte Carlo it is possible to sample
particle paths with correct weight in predefined temperature and collect data
enough for convergence of expectation values of relevant operators.

This approach based on imaginary time propagation is called Path Integral
Monte Carlo method (PIMC). David Ceperley and coworkers have carried out
seminal development work and a sizable number of PIMC simulations of various
many-particle systems ranging from superfluid He Pollock and Ceperley (1984)
and neutron matter Blasdell, Ceperley, and Simmons (1993) to electrons and
hydrogen in extreme conditions PIerleoni  _et al._ (1996), including both
bosonic and fermonic particles.

In recent years, one of the authors of this tutorial (TTR) has significantly
contributed to taking PIMC simulations to new though simple quantum particle
systems, such as small atoms Leino, Nieminen, and Rantala (2006); Leino,
Kylänpää, and Rantala (2007), molecules Kylänpää, Leino, and Rantala (2007);
Kylänpää, Rantala, and Ceperley (2012); Tiihonen, Kylänpää, and Rantala (2015,
2017), a chemical reaction Kylänpää and Rantala (2011) and quantum dots Leino
and Rantala (2007, 2004); Tiihonen  _et al._ (2016), with the ultimate goal of
providing a more accurate description of their electronic structure and
related properties including many-body effects, and how they change with
temperature Leino, Kylänpää, and Rantala (2007); Kylänpää and Rantala (2009,
2010). In this context, PIMC has proven to be a very reliable and excellent
method to calculate the electric polarisabilities of atomic and molecular
systems, leading therefore to an accurate estimation of the optical properties
of both individual small quantum systems and collections of them, in the form
of dilute gases Tiihonen, Kylänpää, and Rantala (2017, 2018, 2016).

A different approach, based on a real, rather than imaginary, time path
integral (RTPI) has been recently proposed as a way to describe the full
quantum dynamics of a quantum system at zero-Kelvin, and to also characterise
the evolution of its eigenstates Ruokosenmäki and Rantala (2015); Leino and
Rantala (2004); Ruokosenmäki  _et al._ (2017); Gholizadehkalkhoran,
Ruokosenmäki, and Rantala (2018); Tiihonen, Kylänpää, and Rantala (2018);
Ruokosenmäki and Rantala (2018). A combination of PIMC and RTPI therefore
gives the possibility to have a comprehensive tool to study the properties of
complex systems and their classical and quantum evolution. This feature in
particular might prove to be very useful not only in chemistry and condensed
matter physics, where this technique fluorished in the past decades, but also
as a viable mean to understand and design the properties of materials of
interest for photonics.

A fully integrable simulation platform, that allows control of both electronic
and photonic properties of matter _exactly_ , without the necessity to revert
to approximations or effective theories, in fact, would constitute a
tremendous resource towards the optimisation of integrated photonic systems.

It is interesting to notice, that throughout the last 30 years, path integrals
have been used to describe several problems in classical and quantum optics,
such as the propagation of light in gradient index media Gómez-Reino and
Liñares (1987), the estimation of the channel capacity of classical and
quantum fiber-based communication networks Reznichenko, Terekhov, and Turitsyn
(2017), parametric amplification Hillery and Zubairy (1984, 1982), light-
matter interaction beyond the rotating wave approximation Zaheer and Zubairy
(1988), decoherence and dephasing in nonlinear spectroscopy Tanimura and
Mukamel (1993), and the effect of retardation in radiative damping Chernyak
and Mukamel (1994). Path integrals have also been employed to link the
nonparaxial propagation of light with different models for quantum gravity
Braidotti  _et al._ (2018). All these works share the common thread of
employing nonrelativistic path integrals to calculate the propagator (_i.e._ ,
the Green’s function) of the electromagnetic field in different contexts, and
use this information to solve the problem at hand. A different approach, based
on path integrals in quantum field theory and Feynman diagrams, has been
recently introduced as a viable way to handle classical Bechler (1999) and
quantum Difallah, Szameit, and Ornigotti (2019) optical phenomena in arbitrary
media.

However, the benefits of path integrals in photonics, namely their ability to
calculate both the properties of matter and its interaction with light in an
exact way, without the need of any approximation both on the matter and light
side, and the new physical insight that this could bring to photonics, remains
to date uncharted territory.

In this tutorial, we aim at introducing the concept and methods of path
integrals to the reader unfamiliar with the field, and to provide researchers
in optics and photonics with a reference point for both analytical and
numerical methods involving path integrals, with the hope that this will
provide a powerful and practical toolkit that could be used in the future to
tackle challenging problems in photonics. An accurate description of the
diverse interactions between light and matter, emerging from the interplay of
fundamental particles and fields, calls naturally for the use of quantum
physics, and its degree of complexity grows very quickly. Path integrals are a
natural way to study these interactions, and they actively take advantage of
the complexity of the problem. Using them in photonics might then lead to
novel methods to exploit complicated light-matter dynamics in photonic
systems.

The tutorial is split into two main parts: Part 1, comprising Sections II
through VI, covers the basics of the path integral approach in physics, and
presents examples on how they can be used to solve problems in classical and
quantum optics, as well as how to employ PIMC to determine optical properties
of materials. In particular, Sect. II covers the fundamentals of the path
integral approach in physics, treating core concepts such as the principle of
least action, classical and quantum probabilities, and a brief description of
the mathematics of integration over infinite number of paths. Basic examples,
comprising the dynamics of a free quantum particle, the quantum harmonic
oscillator, and diffraction from a double slit are covered in Sect. III. To
conclude Part I, two examples of the use of path integrals in classical and
quantum optics are presented, namely the propagation of light in an
inhomogeneous medium and how this could be related to the physics of a
harmonic oscillator with a time-dependent frequency, in Sect. IV, and the
investigation of degenerate parametric down conversion presented in Sect. V
(based on Ref. 37), respectively. Finally, Sect. VI briefly discusses how PIMC
can be used to predict the optical properties of matter, and presents some
perspectives on the use of this computational resource for photonics.

Part 2, on the other hand, including Sections VII through XI, deals with the
basics of path integrals in quantum field theory (QFT), and presents an
application of such framework to the case of the dynamics of the
electromagnetic field in arbitrary media. In particular, Sect. VII briefly
introduces the concept of path integral for quantum fields, and makes use of
the simple case of a scalar field as an example to calculate the relevant
quantities and establish the formalism. After having done that, Sect. VIII
discussed how to include nonlinear interactions in the formalism, and
introduces Feynman diagrams. The results from these two sections are then
intuitively and qualitatively generalised for the case of a vector field in
Sect. IX, as a reference point for Sect. X, where these results are applied to
the particular case of an electromagnetic field propagating in a dispersive
medium of arbitrary shape. The last section of Part II, namely Sect. XI then
presents two explicit examples, of how path integrals can be used in quantum
optics. The first example presents how to describe the onset of spontaneous
parametric down conversion (SPDC) in lossy media through path integrals, while
the second example deals with the calculation of the rate of spontaneous
emission of a quantum emitter surrounded by a dispersive medium.

In the spirit of the educational purpose of a tutorial, and given the
mathematical complexity of path integrals especially concerning the concepts
introduced in Part 2, we also provide, in Appendix A, a step-by-step guide on
how to deal with path integral calculations for the explicit case of the
electromagnetic field in arbitrary media. We hope this would serve as a good
reference and guide to better understand the techniques and methods presented
below.

Finally conclusions and future perspectives are then given in Sect. XII.

##  II Fundamentals of Path Integrals

###  II.1 Probability Amplitudes: Classical vs. Quantum

Figure 1: Pictorial representation of the evolution of a classical (left) and
quantum (right) system from an initial state at time $t_{a}$, to a final state
at time $t_{b}$ over a certain potential landscape. (a) While many different
paths link the initial and the final state of the classical system, once the
initial conditions have been fixed, only the path with minimal action (red
line in the figure) is the actual one undertaken by the system. (b) For a
quantum system, on the other hand, the classical path (blue tube) minimising
the action is interpreted as the path along which the interference of all the
contributing paths linking $t_{a}$ and $t_{b}$ is maximum.

Quantum physics is an abstract theory, whose specific features beyond
classical physics are typically only spectroscopically observable. A good
starting point to find the underlying differences between the two seemingly
different worlds of classical and quantum physics, is represented by their
different interpretation of the concept of probability. This, in fact, turns
out to be a direct manifestation of the wave nature of quantum particles, and
thus, the fundamental issue that we need to incorporate in the study of the
dynamics of quantum systems.

The necessity of a change in viewpoint concerning probability, and the
consequent definition, for quantum physics, of a complex-valued probability
amplitude, emerges very clearly within the context of the least action
principle. Let us consider the situation depicted in Fig. 1, where a classical
[panel (a)] and a quantum [panel (b)] particle are evolving from an initial
time $t_{a}$ to a final time $t_{b}$.

For the classical, deterministic, system, although many different paths
joining $t_{a}$ with $t_{b}$ are available, only the stationary paths with
least action [red line in Fig. 1 (a)] gives a significant contribution to its
dynamics, and determine, ultimately, its equations of motion.

For a quantum system, on the other hand, its intrinsic wave nature (and,
ultimately the uncertainty principle) prevent it for following one single
path, and the classical path [blue tube in Fig. 1 (b)] must be interpreted as
the one with maximum constructive interference coming from all possible paths.

To understand this better, let us consider a classical particle first
propagating from a point $x_{a}$ to a point $x_{c}$, and subsequently to a
third point $x_{b}$. If we denote with $P(c,a)\equiv P(x_{a}\rightarrow
x_{c})$ the probability for the particle to propagate from $x_{a}$ to $x_{c}$,
and, similarly, with $P(b,c)\equiv P(x_{c}\rightarrow x_{b})$ the probability
for the particle to propagate from $x_{c}$ to $x_{b}$, the (classical)
conditional probability for the particle to propagate from $x_{a}$ to $x_{b}$
by going through $x_{c}$ reads

$P(b,a)=\sum_{c}P(b,c)P(c,a),$ (1)

where the summation over $c$ takes into account all possible alternatives for
the intermediate state $c$. Note, that the definition of conditional
probability given above, _i.e._ , $P(b,a)$, differs slightly from the usual
one, which reads $P(a|b)$. This notation, however, is fully equivalent with
the traditional one, and will turn out to be of more practical use for the
purpose of this work.

The above definition can be readily generalised to the case of continuous
variables, _i.e._ , to probability densities, by promoting $P(a)$ and $P(b)$
to probability density functions, and to interpret $a$ and $b$ as two sets of
coordinates $\\{x_{a}\\}$ and $\\{x_{b}\\}$ for the particle to occupy at a
given time $t_{a}$ and $t_{b}$, respectively 222we implicitly assumed that the
two set of coordinates span a one dimensional space, for simplicity. Extension
to higher dimensions can be, in fact, easily carried out.. In this case, then,
the summation over all possible alternatives $c$ in Eq. (1) becomes an
integral over the set of coordinates $\\{x_{a}\\}$, _i.e._ ,

$P(b)=\int_{a}P(b,a)P(a)\,{\rm d}x_{a},$ (2)

where the subscript $a$ on the integral indicates that the integration has to
take into account all the possible values of the integrating coordinate
$x_{a}$.

Now, let us extend the concepts introduced above to the case of a quantum
particle with wave nature. To do that, let us first rewrite Eq. (1) for the
probability amplitude $\psi$ associated to the quantum particle as

$\psi(b,a)\equiv K(b,a)=\sum_{c}\psi(b,c)\psi(c,a).$ (3)

The probability amplitude $\psi(b,a)$ defined above constitutes the basic
quantity from which the dynamics of a quantum particle can be derived, and it
is usually referred to in literature as the kernel (or propagator, or Green’s
function) of the quantum system at-hand. The conventional symbol for that in
path integral language is $K(b,a)$, and we then adopt this notation for the
rest of this manuscript.

The kernel defined above has a simple physical interpretation. In fact, it can
be thought as the impulse response of the system at-hand Feynman and Hibbs
(1965); Byron and Fuller (1992). Moreover, if we know the probability
amplitude of the system at a given initial state $\psi(a)$, we can immediately
evolve it to a final state characterised by a probability amplitude $\psi(b)$
using the following relation

$\psi(b)=\int_{a}K(b,a)\psi(a)\,{\rm d}x_{a}.$ (4)

Finally, the experimentally observable classical probability distributions are
found as squares of the absolute values of the probability amplitudes,
$P(a)=|\psi(a)|^{2}$ and $P(b)=|\psi(b)|^{2}$ at times $t_{a}$ (initial state)
and $t_{b}$ (final state), respectively. With this definition, the probability
amplitudes appearing above can be readily interpreted as the wave function of
the quantum system. Equation (3), in particular, hints at the interpretation
of the wave function of a quantum system as the sum (or, better said,
interference) of all possible paths linking the initial and final states of
the considered evolution.

###  II.2 Lagrangian, Action and Path Integral

Contrary to the canonical formulation of quantum mechanics, which bases its
premises on the Hamiltonian function $H$ of the system, and therefore on the
concept of total energy Sakurai (1994), the path integral formalism starts
from the Lagrangian function, generally defined as the difference between
kinetic and potential energy of the system Arnold (1997), _i.e._ , $L=T-V$.
For this reason, the path integral formalism is often referred to as the
_third formulation_ of quantum mechanics, with the first being the matrix
mechanics developed by Heisenberg, Born and Jordan in 1925 Heisenberg (1925);
Born and Jordan (1925); Jordan (2006), while the second one being the familiar
Hamiltonian formulation developed by Schrödinger in 1926 Schrödinger (1926).

Although the Hamiltonian and Lagrangian formulation of physical problems are
essentially, from the point of view of physical meaning, equivalent, the
latter is more elegant, and, thanks also to its natural appearance in the path
integral formalism, has been adopted as the natural framework for more
complicated theories, such as QFT Brown (1994); Das (2006), particle physics
Kaku (1993), and string theory Green, Schwarz, and Witten (2012).

For simplicity of notations, we proceed in using the one dimensional space
with the coordinate $x$, but the generalisation to three dimensions is
trivial. Then, for a particle with mass $m$ in motion on the path $x(t)$ with
velocity $\dot{x}(t)$ in a potential $V(x,t)$ the classical Lagrangian is

$L(\dot{x},x,t)={1\over 2}m\dot{x}^{2}-V(x,t).$ (5)

On a way to both finding the classical equation of motion and concurrently
incorporating the wave nature of dynamics of the particle, on the path $x(t)$
from $t_{a}$ to $t_{b}$ we define the action

$S(b,a)=\int_{t_{a}}^{t_{b}}\,dt\,L(\dot{x},x,t).$ (6)

In Lagrangian mechanics the action is a parameter related to the path length,
whose optimisation will lead to the equations of motion. This procedure,
_i.e._ , optimising the path $x(t)$ such that $\delta S=0$ is called “the
principle of least action", and leads to the Euler–Lagrange equation Landau
and Lifchitz (1996); Arnold (1997)

$\frac{\rm d}{{\rm d}t}\left(\frac{\partial
L}{\partial\dot{x}}\right)-\frac{\partial L}{\partial x}=0.$ (7)

If we now substitute Eq. (5) into the equation above, we find the following
differential equation

$m\ddot{x}=-\frac{\partial V(x)}{\partial x},$ (8)

which we immediately recognise as Newton’s classical equation of motion. This
constitutes the fundamental ingredient for defining the classical and quantum
probabilities as in Eqs. (1) and (3). However, for the dynamics of a classical
system, the integral approach above is redundant, and the traditional approach
based on Newton’s equation of motion is favourable one in most cases.
Similarly, in quantum mechanics Schrödinger equation is the best approach for
simple problems. However, there are several sophisticated cases where Eqs. (3)
and (4) are more practical and easy to handle. For these more complicated
scenarios, therefore, we need to find the kernel $K(b,a)$, and the right way
to do that is to incorporate in the “classical" path-based approach above, the
information of the wave nature of quantum particles, _i.e._ , to introduce
interference between paths.

Following the ideas of Dirac Dirac (1933) and Feynman Feynman (1948); (2005)
(ed.), we consider the action $S$ as a (classical) measure of the path length,
and Planck’s constant $h$ as the wave length. Then, we can assign a wave
$\exp{\left[i2\pi S(x(t),a)/h\right]}$ to the path $x(t)$ and follow the phase
of the waves to find the interference effects. In particular, we want to take
into account the contributions of all possible waves of the form
$\exp{\left[iS(b,a)/\hbar\right]}$ to the probability amplitude $K(b,a)$.

The sum (or integral) over the contributions of all possible paths is called
the path integral, _i.e._ ,

$K(b,a)=\int_{a}^{b}\,Dx(t)\,\exp{\left[\frac{i}{\hbar}S(b,a)\right]},$ (9)

where the notation ${Dx(t)}$ indicates integration over all paths from
$a=(x_{a},t_{a})$ to $b=(x_{b},t_{b})$, which, following Ref. 2, can be
defined as

$\int_{a}^{b}\,Dx(t)=\lim_{\varepsilon\rightarrow
0}\frac{1}{A(\varepsilon)}\int\,\frac{dx_{1}}{A(\varepsilon)}\,\int\,\frac{dx_{2}}{A(\varepsilon)}\,\cdots\,\int\,\frac{dx_{N-1}}{A(\varepsilon)},$
(10)

where $A(\varepsilon)$ is a suitable normalisation factor that ensures the
limit to properly converge (an example of it is given in the next section, but
the explicit for of $A(\varepsilon)$ might change depending on the problem at-
hand), and $\varepsilon$ comes from the discretisation of the time interval in
the action $S(b,a)$ in $N$ finite points, i.e., $t_{b}-t_{a}=N\varepsilon$,
with $t_{0}=t_{a}$, and $t_{N}=t_{b}$. This, on the other hand, implies a
discretisation on the paths $x(t)$, which now are defined as
$x_{k}=x(t_{k})=x(t_{0}+k\varepsilon)$, with $x_{0}=x_{a}=x(t_{a})$ and
$x_{N}=x_{b}=x(t_{b})$. This discretisation procedure allows us to consider
each of the $N-1$ integrals above as standard Riemann/Lebesgue integrals in
the variable $dx_{k}$. Then, when all $N-1$ integrals have been computed, the
limit $\varepsilon\rightarrow 0$, together with the correct definition of the
normalisation factor $A(\varepsilon)$, ensures convergence of the path
integral $\int_{a}^{b}\,Dx(t)$, and justifies its definition as integration
over all possible paths. Notice, moreover, that the same line of reasoning
will allow us, in Sect. VII to define the path integral for fields.

The explicit form of the path integral above is usually defined by the form of
the potential term $V(x)$ appearing in the Lagrangian (5). For some potential
functions there are analytical, closed-form, exact solutions to the Eq. (9),
but one needs to prepare for numerical methods with possible approximations in
more general cases, such as multi-dimensional or many-body problems. In what
follows, we present some of the exact propagators.

The above result, obtained for a simple one dimensional system, can be readily
extended to three dimensional and many particle systems. While the former is
straightforward to work out, in the latter case, the quantum statistics of
fermions and bosons needs to be taken into account explicitly, which makes the
problem of finding the correct generalisation of the path integral to the many
body case less trivial Fetter and Walecka (2003).

In addition to that, the path integral approach also allows for an easy way to
simulate the evolution of the density matrix for finite temperature
equilibrium systems. An example of this will be discussed below in Sect. V.

##  III Basic Examples for Quantum Particles

For the evaluation of the wave function from the integral equation (4) we need
to calculate explicitly the kernel from the path integral, _i.e._ , Eq. (9).
In this section we therefore consider the simplest kernels, with examples
following the book of Feynman and Hibbs Feynman and Hibbs (1965).

In case the integrand is an exponential of a quadratic function (Gaussian
integral) the kernel can always be evaluated using recursively the basic
Gaussian integral formula Byron and Fuller (1992)

$\int\,dx\,e^{-ax^{2}+bx}=\sqrt{\frac{\pi}{a}}\,e^{b^{2}/4a}.$ (11)

Another useful result to keep in mind is the fact, that for given $x_{a}$ and
$x_{b}$ endpoints, the contributions from other than the classical path
interfere destructively and vanish. Thus, only the classical action
$S_{cl}(x_{b},x_{a})$ contributes, leading to the major simplification for the
kernel, _i.e._ ,

$K(x_{b},x_{a})=F(t_{b}-t_{a})\exp{\left[\frac{i}{\hbar}S_{cl}(x_{b},x_{a})\right]}.$
(12)

As long as the action involves path variables only up to the second order,
therefore, the exact propagators in the form above can be factored out from
the path integral, leaving at most to calculate a prefactor of the form
$F(t_{b}-t_{a})$.

To illustrate how one arrives at the result above, let us consider a general
quadratic Lagrangian in $x$ and $\dot{x}$ of the form

$L(x,\dot{x},t)=a\dot{x}^{2}+b\dot{x}x+cx^{2}+d\dot{x}+ex+f,$ (13)

where $\\{a,b,c,d,e,f\\}$ are (arbitrary, but well-behaved) time-dependent
coefficients. The action corresponding to this Lagrangian is then the integral
of Eq. (13) with respect to time between two fixed end points $t_{a}$ and
$t_{b}$, as given by Eq. (6).

Let us now assume, that $x_{cl}(t)$ is the classical path between the
specified end points, i.e., the path for which $\delta S=0$ holds. We can then
represent $x(t)$ in terms of deviations from the classical path $x_{cl}(t)$ by
introducing the function $y(t)$ as

$x(t)=x_{cl}(t)+y(t).$ (14)

This substitution means, that instead of defining a points on the path by its
distance from and arbitrary coordinate axis, we measure instead the deviation
$y(t)$ from the classical path. Moreover, at each time $t\in[t_{a},t_{b}]$,
the variables $x$ and $y$ differ only by the constant $x_{cl}(t)$, and
therefore $dx_{i}$ = $dy_{i}$ for each point $t_{i}$. In general, then, it
follows that $Dx(t)=Dy(t)$. Notice, that as a consequence of Eq. (14),
$y(t_{a})=0=y(t_{b})$, as at the endpoints the path $x(t)$ coincides with the
classical path $x_{cl}(t)$.

If we then use the change of variables defined in Eq. (14), the Lagrangian
(13) can be written as the sum of three terms as follows

$L(x,\dot{x},t)=L_{cl}+L_{y}+L_{mix},$ (15)

where $L_{cl}$ ($L_{y}$) is just Eq. (13) with
$\\{x,\dot{x}\\}\rightarrow\\{x_{cl},\dot{x}_{cl}\\}$
($\\{x,\dot{x}\\}\rightarrow\\{y,\dot{y}\\}$), and

$L_{mix}=\left(2a\dot{x}_{cl}+bx_{cl}\right)\dot{y}+\left(b\dot{x}_{cl}+2cx_{cl}\right)y.$
(16)

Similarly, the action can be then written as the sum of three terms, namely
the classical action $S_{cl}(b,a)$, the action relative to the deviation
$y(t)$, i.e., $S_{y}(b,a)$, and the mixed action $S_{mix}(b,a)$. Because of
the fact that $y(t_{a})=0=y(t_{b})$, however, all the terms which contain
linear terms in $y$ result in a vanishing integral. Thus, only the second-
order terms in $y$ give rise to a nonzero contribution to the total action,
which can now be written as

$S(b,a)=S_{cl}(b,a)+\int_{t_{a}}^{t_{b}}\left(a\,\dot{y}^{2}+b\,\dot{y}y+c\,y^{2}\right)dt.$
(17)

Notice how $S_{cl}(b,a)$ does not depend on the deviation $y(t)$, and
therefore the corresponding exponential can be treated as a constant, with
respect to the path integration $Dy(t)$. The Kernel can then be written in the
following form

$\displaystyle K(b,a)$ $\displaystyle=$
$\displaystyle\int_{0}^{0}\,Dy(t)\,\exp{\left[\left\\{\frac{i}{\hbar}\int_{t_{a}}^{t_{b}}(a\,\dot{y}^{2}+b\,\dot{y}y+c\,y^{2})dt\right\\}\,\right]}$
(18) $\displaystyle\times$
$\displaystyle\exp{\left[\frac{i}{\hbar}S_{cl}(b,a)\right]},$

where the notation $\int_{0}^{0}\,Dy(t)$ is reminiscent of the fact that all
the paths $y(t)$ obey the boundary condition $y(t_{a})=0=y(t_{b})$. The path
integral above, then, can be written as a function of the time interval
$(t_{b}-t_{a})$ solely, i.e.,

$F(t_{b}-t_{a})=\int_{0}^{0}\,Dy(t)\,\exp{\left[\left\\{\frac{i}{\hbar}\int_{t_{a}}^{t_{b}}(a\,\dot{y}^{2}+b\,\dot{y}y+c\,y^{2})dt\right\\}\,\right]}.$
(19)

This, ultimately, allows us to write the kernel in the following simplified
form

$K(b,a)=F(t_{b}-t_{a})\exp{\left[\frac{i}{\hbar}S_{cl}(x_{b},x_{a})\right]},$
(20)

which is equivalent to that of Eq. (12).

###  III.1 Path Integrals for a Free Quantum Particle

Figure 2: Real (blue, solid line) and imaginary (red, dashed line) parts of
the kernel $K(x_{b},x_{a})$ of a free quantum particle, as given by Eq. (30),
as a function of $x=x_{b}-x_{a}$, for a fixed time interval $t_{b}-t_{a}=T$.
For this plot, we assumed $m=1=\hbar$, which corresponds to measure time in
square meters, rather than seconds. Moreover, $T=1$ has also been used.

The first example concerns the simplest quantum system, _i.e._ , a quantum
particle of mass $m$, freely propagating without experiencing any interaction.
Following the assumptions made in the previous section, we discuss the case of
a one dimensional free particle. The generalisation to an arbitrary number of
dimensions can be readily done, since the dynamics of a free quantum particle
in $D$ dimensions can be seen as the product of the independent evolution of
$D$ one dimensional particles Dirac (1947).

The Lagrangian of a free particle of mass $m$ is given by

$L(x,\dot{x})=\frac{m\,\dot{x}^{2}}{2},$ (21)

and the equation of motion deriving from the Euler-Lagrange equations (7) is
simply $m\,\ddot{x}=0$.The correspondent action
$S=\int_{t_{a}}^{t_{b}}L(\dot{x},x,t)dt$ can be readily calculated explicitly
by means of part integration, and has the following form

$S=\frac{m(x_{b}-x_{a})^{2}}{2(t_{b}-t_{a})},$ (22)

where $x_{a,b}=x(t_{a,b})$. To calculate the path integral for a free
particle, we need to consider all the possible paths the particle takes from
the initial state $(t_{a},x_{a})$ to the final state $(t_{b},x_{b})$. To do
that, we simply divide the time interval $T=t_{b}-t_{a}$ into $N$ smaller
intervals of length $\varepsilon=t_{i+1}-t_{i}$ (such that $T=N\varepsilon$),
calculate the action $S_{i}$ corresponding to the particle evolution within
each single interval, such that

$S=\sum_{i=1}^{N}S_{i}=\sum_{i=1}^{N}\frac{m}{2\varepsilon}(x_{i}-x_{i-1})^{2},$
(23)

with $x_{0}\equiv x_{a}$ and $x_{N}\equiv x_{b}$. We then substitute this
result into Eq. (9) and evaluate the path integral over the set of $N$
distinct trajectories, i.e,

$\int_{a}^{b}\,Dx(t)\hskip 14.22636pt\rightarrow\hskip
14.22636pt\frac{1}{A(\varepsilon)^{N}}\int\,dx_{1}\int\,dx_{2}\,\cdots\,\int\,dx_{N-1},$
(24)

where $A(\varepsilon)=\sqrt{(2\pi i\hbar\varepsilon)/m}$ is a factor included
to ensure the integral to converge. This factor, however, is not merely a
normalisation factor, since it is complex, and therefore it contributes to the
overall phase of the path integral. This is discussed in great detail in Ref.
2. Finally, we take the limit $N\rightarrow\infty$, to arrive at the following
expression for the propagator of a free particle

$K(x_{2},x_{1})=\lim_{N\rightarrow\infty}\frac{1}{A(\varepsilon)^{N}}\int\,dx_{1}\cdots\int\,dx_{N}\exp{\left(\frac{i}{\hbar}\sum_{i=1}^{N}S_{i}\right)}$
(25)

The integrals appearing above are Gaussian in the variables $x_{i}$, and can
be then readily calculated one after another. To see how, let us first
calculate explicitly the integral with respect to $x_{1}$. Once we have this
result, we can perform the other integrations in cascade in the same manner.

First, notice that the relevant term in the action, depending explicitly on
$x_{1}$ (_i.e._ , those obtained by setting $i=1$ and $i=2$ in the expression
above) give rise to the following term

$\displaystyle\frac{i}{\hbar}S_{1}$ $\displaystyle=$
$\displaystyle\frac{im}{2\hbar\varepsilon}\left[(x_{2}-x_{1})^{2}+(x_{1}-x_{0})^{2}\right]=\frac{im}{2\hbar\varepsilon}(x_{2}^{2}+x_{0}^{2})$
(26) $\displaystyle+$
$\displaystyle\frac{im}{\hbar\varepsilon}\left[x_{1}^{2}-(x_{2}+x_{0})x_{1}\right],$

and that, in particular, the first term does not depend on the integration
variable $x_{1}$. Integrating the above quantity with respect to $x_{1}$ then
gives

$\displaystyle\frac{m}{2\pi
i\hbar\varepsilon}\int\,dx_{1}\,\exp{\left(\frac{i}{\hbar}S_{1}\right)}=\exp{\left[\frac{im}{2\hbar\varepsilon}(x_{2}^{2}+x_{0}^{2})\right]}$
(27) $\displaystyle\times$ $\displaystyle\frac{m}{2\pi
i\hbar\varepsilon}\int\,dx_{1}\,\exp{\left\\{\frac{im}{\hbar\varepsilon}\left[x_{1}^{2}-(x_{2}+x_{0})x_{1}\right]\right\\}}$
$\displaystyle=$ $\displaystyle\sqrt{\frac{m}{2\pi
i\hbar(2\varepsilon)}}\exp{\left[\frac{im}{2\hbar(2\varepsilon)}(x_{2}-x_{0})^{2}\right]},$

where to pass from the second to the third line we have employed the change of
variables $X=ix_{1}$ and used Eq. (11).

Next, we take into account the terms in the action depending explicitly on
$x_{2}$, _i.e._ , $S_{2}=S_{1}+(x_{3}-x_{2})^{2}$, and we integrate with
respect to $x_{2}$. We can do so, by simply taking the result of the integral
of $S_{1}$ given in Eq. (27) and multiply it by

$\sqrt{\frac{m}{2\pi
i\hbar\varepsilon}}\exp\left[\frac{im}{2\hbar\varepsilon}(x_{3}-x_{2})^{2}\right],$
(28)

and integrate again, this time over $x_{2}$. The result is similar to that of
Eq. (27), except that $(x_{2}-x_{0})^{2}$ becomes $(x_{3}-x_{0})^{2}$, and
$2\varepsilon\rightarrow 3\varepsilon$. It is now clear, that we can solve the
set of $N$ integrals in Eq. (25) by recursively applying terms of the form
(28) and then performing Gaussian integration with respect to the variable
$x_{i}$. After $N-1$ steps, we are left with the following result

$\sqrt{\frac{m}{2\pi
i\hbar(N\varepsilon)}}\exp\left[\frac{im}{2\hbar(N\varepsilon)}(x_{N}-x_{0})^{2}\right].$
(29)

If we now notice, that $N\epsilon=T=t_{b}-t_{a}$, it is easy to see that the
limit operation in Eq. (25) can be readily performed, leading us to the final
form for the path integral of a free particle, _i.e._ ,

$K(x_{b},x_{a})=\sqrt{\frac{m}{2\pi
i\hbar(t_{b}-t_{a})}}\exp{\left[\frac{im(x_{b}-x_{a})^{2}}{2\hbar(t_{b}-t_{a})}\right]}.$
(30)

The functional form of the real and imaginary parts of the kernel above for a
constant time interval $t_{b}-t_{a}=T$ are shown in Fig. 2. It is worth to
comment this result a bit. From it, in fact, we see how the quantum dynamics
of a free particle (but, more in general, of an arbitrary quantum system)
couples with its classical dynamics described by the action (22). However,
this result also allows to shed light on the essential difference between
quantum and classical dynamics. While the classical principle of least action
localised the propagation of a particle on a specific trajectory $x(t)$
(_i.e._ , the path of minimal action), the propagator $K(x(t),x_{a})$
describing the propagation of a particle from the initial state $x_{a}$ along
the classical trajectory $x(t)$, yields instead a complex wave function
$\psi(x,t)$ delocalising in time during propagation (obviously within the
constraints of the Heisenberg principle).

###  III.2 Refraction of Photons at an Interface

Figure 3: Schematic representation of a ray of light being refracted from a
planar interface. A light source located at $A$ emits a photon, which gets
detected at point $B$ by a suitable detector. Between $A$ and $B$, a planar
interface $C$ is placed, separating the two different media in which the
photon propagates. In the language of path integrals, refraction can be
interpreted as the necessary change in velocity of the photon along the path
where the action has the minimal value, with respect to the contact point
$x_{c}$ on the planar interface $C$. Figure 4: Left panel: pictorial
representation of the double-slit experiment. A quantum particle is emitted
from a source located in $A$ and reaches a detector (or a distribution
thereof) located at point $B$ by passing through a screen $C$, on which two
slits, $C_{1,2}$ have been carved. Right panel: sketch of the various possible
probability distributions that could be observed at the detector plane $B$.
(a) Interference pattern $P$ revealing the wave nature of the quantum
particle. (b) Probability distribution $P_{2}$ obtained by closing the slit
$C_{1}$. (c) Probability distribution $P_{1}$ obtained by closing the slit
$C_{2}$. (d) Estimation of the total probability distribution $P_{1}+P_{2}$,
obtained by simply summing the distributions obtained in the cases (b) and
(c), as would be valid for a classical particle, where no interference occurs.

With the above notations, consider now the space divided in two parts by a
planar interface $C$, and let $A$ and $B$ be two points located at opposite
sides of the interface, as it is shown in Fig. 3. Then, we assume constant but
different potentials at opposite sides of the interface. At both sides we
expect the path of the photons to be those of a free particle. With this
assumption, we can project the path connecting the points $A$ and $B$ onto a
one dimensional subspace, for simplicity. Notice, however, that despite we
will perform the calculations in this one dimensional subspace, we still need
to think in terms of three dimensional space when considering deflection of
light rays from the interface $C$.

Without loss of generality, we can assume that total energy is conserved in
both separate sides of the interface, and while passing through it. However,
since the two sides might have different values of the potentials (_i.e._ ,
different refractive indices), the particle passing through the interface $C$
needs to change its velocity from $\dot{x}_{a}$ to $\dot{x}_{b}$, to account
for this variation in potential energy. The change of velocity, moreover,
occurs at some position $x_{c}$ on the interface.

As shown in Fig. 3, in case the straight line from $A$ to $B$ is not
perpendicular to the interface, the observed path becomes deflected at
$x_{c}$.

Now, we can write the action for the path as

$\displaystyle S_{cl}(x_{b},x_{a})$ $\displaystyle=$ $\displaystyle
S_{cl}(x_{b},x_{c})+S_{cl}(x_{c},x_{a})$ (31) $\displaystyle=$
$\displaystyle\frac{m}{2}\left[\dot{x}_{b}\left(x_{b}-x_{c}\right)+\dot{x}_{a}\left(x_{c}-x_{a}\right)\right],$

where the coordinate $x_{c}$ is a free parameter to optimise following the
principle of least action. For the case of a photon refracting from interface
$C$, therefore, the "minimum optical path length" from $A$ to $B$ is found by
optimising the coordinate $x_{c}$. Following Eq. (12), the kernel for the
refraction of a quantum particle from a planar interface $C$ is given, up to
an inessential constant $F\equiv F(t_{b}-t_{a})$, by

$K(x_{b},x_{a})=F\exp{\left\\{\frac{i\,m}{2\hbar}\left[\dot{x}_{b}\left(x_{b}-x_{c}\right)+\dot{x}_{a}\left(x_{c}-x_{a}\right)\right]\right\\}}.$
(32)

We leave as an exercise for the reader to figure out the explicit expression
of the normalisation constant $F\equiv F(t_{b}-t_{a})$.

In the equations above, the two different constant velocities $\dot{x}_{a}$ to
$\dot{x}_{b}$ of the photon follow from two different refractive indices. In
optics, the optimisation of the optical path length is called Fermat’s
principle. If we consider the problem from a geometrical optics perspective,
in fact, it is easy to see how the optimisation of $x_{c}$ leads directly to
the celebrated Snell’s law of refraction Luneburg (1965), _i.e._ ,

$n_{a}\sin\theta_{a}=n_{b}\sin\theta_{b},$ (33)

where $n_{a,b}=1/\dot{x}_{a,b}$ are the refractive indices of the two media
separated by interface $C$ (that can be expressed as the inverse velocity of
the particle in each side of $C$), and $\theta_{a,b}$ the angle the optical
rays emerging from $A$ and $B$, resepctively, make with it. The angle
difference between the two sides of the interface $C$ is reminiscent of the
different wavelength that a photon with velocity $\dot{x}_{a}$ and one with
velocity $\dot{x}_{b}$ experiences.

In this simple example, the classical and quantum approaches give identical
explanations for observations. So, we see that the "quantum corrections",
though not absent may give not only small, but even vanishing contribution.

###  III.3 Diffraction from a Double Slit

This is the well known “classic experiment" for demonstrating the wave nature
of light, or, if conducted with quantum particles like electrons, to expose
the wave nature of particle dynamics Feynman (1990). The results of this
section can be then thought as valid for both a photon, or a quantum particle.
Experimentally, to prove the wave nature of light one would need to perform
this experiment with monochromatic light, while mono-energetic particles are
needed to unravel the wave nature of particle dynamics. From a path integral
perspective, the interference pattern typical of such experiments naturally
arises when the possible paths the particle can take to traverse the double
slit are explicitly taken into account.

Our first task is then to construct the kernel for this problem. To this aim,
let us consider the situation depicted in Fig. 4, where a particle is emitted
from a suitable source located at point $A$ and propagates to a detector (or a
collection thereof) placed at $B$, through a screen $C$ with two slits $C_{1}$
and $C_{2}$ carved in it. Without loss of generality, we can assume the
evolution of the particle is that of a free particle, and that the only
potential it encounters is represented by the double-slit structure (which, as
a matter of fact, acts as a transfer function for the particle).

To reach the detector at $B$, the particle takes a time $\tau$ to travel from
$A$ to the screen $C$ (and, in particular, to one of the slits), pass through
the screen, and then arrive at $B$ after a time $T-\tau$, where
$T=t_{b}-t_{a}$ is the total time the particle takes to go from $A$ to $B$.
However, since we don’t know exactly the time at which the particle arrives on
the screen, we need to integrate over all possible times. In other terms,
since we don’t know the exact path the particle will take to reach $C$ from
$A$, we need to integrate over all possible paths it might undertake. With
this in mind, and by remembering that $a=(x_{a},t_{a})$, $b=(x_{b},t_{b})$,
and $c_{1,2}=(x_{c_{1,2}},t_{c_{1,2}})$ represent the position and time at
points $A$, $B$, and $C_{1,2}$, respectively, the kernel for the propagation
of a particle through a two-slit screen is given by

$\displaystyle K(b,a)$ $\displaystyle=$
$\displaystyle\int\,d\tau\,\Big{[}K(c_{1},b;\,\tau)K(a,c_{1};\,T-\tau)$ (34)
$\displaystyle+$ $\displaystyle K(c_{2},b;\,\tau)K(a,c_{2};\,T-\tau)\Big{]},$

where the first term accounts for the particle passing through the slit
$C_{1}$, while the second one for it passing through the slit $C_{2}$. From
the expression above, it is clear that interference must occur, since the
probability to detect the particle at point $B$ is then proportional to
$|K(b,a)|^{2}$.

In fact, we can arrive at the same conclusion by considering the case above in
terms of the probability amplitudes and observed probabilities in Sect. II.1.
In particular, let us have a look at the probability distribution of the
particle observed on the screen at $B$ for different cases, namely the
probability distribution $P$ when interference pattern is observed (both slits
are open) [panel (a)], the probability distribution $P_{2}$ observed when the
slit $C_{1}$ has been closed [panel (b)], the probability $P_{1}$ observed
when the slit $C_{2}$ has been closed [panel (c)], and finally the probability
distribution $P_{1}+P_{2}$ obtained by summing the results of observations in
(b) and (c).

As it can be clearly seen, the probabilities do not sum up, as
$P\not=P_{1}+P_{2}$. If, on the other hand, we first sum the probability
amplitudes, _i.e._ , $\psi=\psi_{1}+\psi_{2}$, and then calculate the
probability distribution as $P=|\psi|^{2}$ we get

$P=|\psi|^{2}=|\psi_{1}+\psi_{2}|^{2}=P_{1}+P_{2}+2\operatorname{Re}\\{\psi_{1}\psi_{2}^{*}\\},$
(35)

where the third term is responsible for adding interference on top of the sum
of probabilities $P_{1}+P_{2}$ in Fig. 4 (d), thus leading to the correct
result of Fig. 4 (a). The reader familiar with wave theory would immediately
recognise that the probability distributions shown in Fig. 4 (a) and (d) occur
as well for classical waves. Specifically, the distribution (a) occurs for the
interference of correlated (i.e., first-order coherent) classical waves, while
distribution (d) results from incoherent classical wave mixing. This is yet
another indicator of the wave nature of quantum particles Feynman (2010).

###  III.4 Path Integral for the Harmonic Oscillator

Figure 5: Real (blue, solid line) and imaginary (red, dashed line) parts of
the kernel $K(x,0)$ of a harmonic oscillator, as defined by Eq. (52), assuming
that the initial condition for the oscillator is $x_{a}=0$, and for a given
time interval $T=1$. Analogously to the convention adopted in Fig. 2,
$m=1=\hbar$ have been assumed to plot this figure as well.

As our next example, we consider a simple harmonic oscillator, described by
the following Lagrangian

$L(x,\dot{x})=\frac{m\dot{x}^{2}}{2}-\frac{m\omega^{2}x^{2}}{2},$ (36)

where $m$ is the mass of the oscillator, and $\omega$ its characteristic
resonance frequency. Using the Euler—Lagrange Eq. (7), the equation of motion
reads $\ddot{x}+\omega^{2}x=0$, and the action
$S=\int_{t_{a}}^{t_{b}}L(\dot{x},x,t)dt$ is explicitly given by Feynman and
Hibbs (1965)

$S=\frac{m}{2T\text{sinc}(\omega T)}\left[(x_{a}^{2}-x_{b}^{2})\,\cos(\omega
T)-2x_{a}x_{b})\right],$ (37)

where $T=t_{b}-t_{a}$, and $\text{sinc}(x)=\sin(x)/x$.

To calculate the path integral for the harmonic oscillator, we take a slightly
different approach, than the one taken above, which allows us to perform
calculations in an easier and more intuitive manner. Let us then assume that
$x_{cl}(t)$ and $\dot{x}_{cl}(t)$ represent the classical path of the
oscillator, and its velocity, respectively. We can then express any other
possible path taken by the oscillator, as a deviation from the classical path,
_i.e._ , $x(t)=x_{cl}(t)+y(t)$ and $\dot{x}(t)=\dot{x}_{cl}(t)+\dot{y}(t)$,
and choose the appropriate boundary conditions on $y(t)$ and $\dot{y}(t)$,
_i.e._ , $y(t_{a})=0=y(t_{b})$ and $\dot{y}(t_{a})=0=\dot{y}(t_{b})$ as
required by the principle of least action.

If we apply this change of variables to the Lagrangian above, we get three
different terms, namely

$L(x_{cl},\dot{x}_{cl})+L(y,\dot{y})+m(\dot{x}_{cl}\dot{y}+\omega x_{cl}y),$
(38)

where the first term has the form given by Eq. (36) and it represents the
Lagrangian for the classical trajectory, the second term is the Lagrangian
(36) for the deviation, and the third term amounts to a vanishing term (for
the Euler—Lagrange equations) and can therefore be neglected Feynman and Hibbs
(1965). Using this result, we can then factor the propagator $K(x_{2},x_{1})$
for the harmonic oscillators into two terms as follows

$K(x_{b},x_{a})=F(T)\exp{\left(\frac{i}{\hbar}S_{cl}\right)},$ (39)

where $S_{cl}=\int_{t_{a}}^{t_{b}}\,dt\,L(x_{cl},\dot{x}_{cl})$, and the
explicit expression of the function $F(T)$, which depends upon the time
interval solely, is

$F(T)=\int\,Dy(t)\,\exp{\left[\frac{i}{\hbar}\int_{0}^{T}\,dt\,\frac{m}{2}(\dot{y}^{2}-\omega^{2}y^{2})\right]}.$
(40)

Because of our choice of boundary conditions for the deviations $y(t)$ and
$\dot{y}(t)$, namely that they must both be zero at the endpoints, we can
significantly simplify the calculation of the path integral above if we allow
the various paths $y(t)$ to be represented as a Fourier series as

$y(t)=\sum_{n=1}^{\infty}a_{n}\sin\left(\frac{n\pi t}{T}\right).$ (41)

The requirement that $y(t_{a})=0=y(t_{b})$ and
$\dot{y}(t_{a})=0=\dot{y}(t_{b})$, in fact, corresponds to say that both
$y(t)$ and $\dot{y}(t)$ are $T$-periodic functions, that can then be
represented in a Fourier series.

This representation of $y$ gives possibility to specify a path through the
coefficients $a_{n}$ instead of values of $y$ at any particular time $t$. This
can be seen as a linear transformation of coordinates, whose Jacobian $J$ is a
dimensionless constant independent of $\omega$, $m$ or $\hbar$. Moreover, we
don’t really need to evaluate $J$ explicitly, since we can always recover the
correct normalisation factor at the end of our calculation, by requiring that

$\lim_{\omega\rightarrow 0}F(T)=\sqrt{\frac{m}{2\pi i\hbar T}},$ (42)

_i.e._ , that in the limit of $\omega\rightarrow 0$, where the Lagrangian of
the harmonic oscillator reduces to that of a free particle, we find the
appropriate normalisation coefficient for a free particle. For this reason, we
omit $J$ from the following calculations, and we restore the correct
normalisation factor only at the very end of them.

Before substituting Eq. (41) into Eq. (40), a couple more assumptions are
needed, in order to easily compute the path integral. First, we truncate the
Fourier series to a finite number $N$, so that

$\int\,Dy(t)\hskip 14.22636pt\rightarrow\hskip
14.22636pt\frac{1}{A(\varepsilon)^{N}}\int\,da_{1}\,\int\,da_{2}\,\cdots\,\int\,da_{N},$
(43)

and we can calculate it recursively as we did for the case of a free particle.
We will then take the limit $N\rightarrow\infty$ at the end of the
calculation.

Putting everything together, we obtain the following expression for the term
$F(T)$, after performing the trigonometric integrals Byron and Fuller (1992),
as

$\displaystyle F(T)$ $\displaystyle=$
$\displaystyle\frac{1}{A(\varepsilon)^{N}}\int\,da_{1}\,\cdots\,\int\,da_{N}\exp\Bigg{\\{}\frac{im}{2\hbar}\frac{T}{2}$
(45) $\displaystyle\sum_{n=1}^{N}\left[\left(\frac{\pi
n}{T}\right)^{2}-\omega^{2}\right]a_{n}^{2}\Bigg{\\}}.$

Notice, once more, that the integrals above are all Gaussian in the
integration variables $a_{n}$. As we did for the case of the free particle,
therefore, we can perform them individually, and then obtain the final result
recursively. The result of a single integration over $a_{n}$ is then given by

$\int\,\frac{da_{n}}{A}\,\exp{\left[\frac{imT}{4\hbar}f_{n}(\omega,T)a_{n}^{2}\right]}=\sqrt{\frac{2}{Tf_{n}(\omega,T)}},$
(46)

where

$f_{n}(\omega,T)=\frac{n^{2}\pi^{2}}{T^{2}}-\omega^{2}.$ (47)

Since there are no linear terms in $a_{n}$ in any of the integrals above, the
final result of the path integration will be proportional to simply the
product of $N$ independent terms $f_{n}(\omega,T)$, one for each value of $n$.
This allows us to write

$F(T)\propto\prod_{n=1}^{N}\frac{1}{\sqrt{f_{n}(\omega,T)}}.$ (48)

We now need to take the limit of the expression above for
$N\rightarrow\infty$. To do that, let us first rewrite the product above in
the following way, using Eq. (47)

$\displaystyle\prod_{n=1}^{N}\frac{1}{\sqrt{f_{n}(\omega,T)}}$
$\displaystyle=$
$\displaystyle\prod_{n=1}^{N}\left(\frac{n^{2}\pi^{2}}{T^{2}}-\omega^{2}\right)^{-1/2}=\prod_{n=1}^{N}\left(\frac{n^{2}\pi^{2}}{T^{2}}\right)^{-1/2}$
(49) $\displaystyle\times$
$\displaystyle\prod_{n=1}^{N}\left(1-\frac{\omega^{2}T^{2}}{n^{2}\pi^{2}}\right)^{-1/2}.$

The first product does not depend on $\omega$ and, together with the Jacobian
and the terms $\sqrt{2/\varepsilon T}$ deriving from the various integrations,
can be collected into an overall normalisation factor $C$. The second product,
on the other hand, admits the following limit

$\lim_{N\rightarrow\infty}\prod_{n=1}^{N}\left(1-\frac{w^{2}T^{2}}{n^{2}\pi^{2}}\right)^{-1/2}=\sqrt{\frac{1}{\text{sinc}(\omega
T)}}.$ (50)

Putting everything together and evaluating $C=\sqrt{m/2\pi i\hbar T}$ from the
free-particle-limit (42) we get the final expression of the term $F(T)$ as

$F(T)=\sqrt{\frac{m}{2\pi i\hbar T\text{sinc}(\omega T)}},$ (51)

and, after substituting this result into Eq. (39), we obtain the final form of
the path integral for the harmonic oscillator to be

$K(x_{b},x_{a})=\sqrt{\frac{m}{2\pi i\hbar T\text{sinc}(\omega
T)}}\exp{\left(\frac{i}{\hbar}S_{cl}\right)},$ (52)

where $S_{cl}$ is given by Eq. (37). The real and imaginary pary of the kernel
$K(b,a)$ for the harmonic oscillator are shown in Fig. 5.

##  IV An Example from Classical Optics: Path Integral Description of Light
Dynamics in an Inhomogeneous Medium

We now take a look at how path integrals can be used to solve problems outside
quantum mechanics, and apply this formalism to describe the propagation of the
electromagnetic field inside a weakly inhomogeneous medium, using the case of
gradient-index (GRIN) media as explicit reference. This problem has been
solved by C. Gõmez-Reino and J. Liñares in 1987 Gómez-Reino and Liñares
(1987). In their work, Gõmez-Reino and Liñares first represent the
electromagnetic field in a GRIN medium as a superposition of optical rays, and
use this picture to calculate the propagator as a path integral over the rays’
trajectories. They then provide an explicit expression for it, parametrised in
the so-called paraxial and field rays of an arbitrary (paraxial) optical
system Luneburg (1965).

Here, we take a different approach, with which we want to show how the free
propagation of light in a medium can be seen as, essentially, the evolution of
a massive quantum particle in a harmonic oscillator potential with a suitably
defined frequency, which, in general, can be $z$-dependent. We will identify
such massive particle with a photon propagating inside the medium, and all the
possible trajectories the particle can take as the possible optical rays
linking the initial ($z=0$) and final ($z=z$) propagation plane in the medium.
We will then calculate the diffraction kernel by means of path integrals,
essentially following the results of Sect. III.4 and show how our calculations
naturally suggest a representation of the diffraction kernel in terms of
Hermite—Gaussian functions.

Without loss of generality, and for the sake of simplicity of exposition, we
consider light propagating in a $1+1$-dimensional GRIN medium, characterised
by the following refractive index profile

$n^{2}(x,z)=n_{0}^{2}[1-g^{2}(z)x^{2}].$ (53)

where $n_{0}=n(0,z)$ is the background index and $g(z)$ is a smooth-enough
function that describes the evolution of the refractive index along the $z$
axis.

###  IV.1 Diffraction Kernel as a Path Integral

Let us first assume paraxial propagation of light in a medium described by the
refractive index given by Eq. (53). In general, if we know the field
distribution at an initial plane $z=0$ to be $E(x_{a},y_{a})$, we can
calculate the field distribution at a plane $z>0$ by means of the diffraction
integral Jackson (1986)

$E(x_{b},z)=\int\,dx_{a}\,K(x_{b},x_{a},z)E(x_{a},0),$ (54)

where $K(x_{b},x_{a},z)$ is the diffraction kernel, _i.e._ , the Green’s
function of the paraxial equation

$2\,i\,k\,n_{0}\frac{\partial K}{\partial
z}=-\left[\frac{\partial^{2}}{\partial x^{2}}+2k^{2}\,n_{0}\,n(x,z)\right]K$
(55)

with the boundary condition that $K(x,x_{a},z)\rightarrow\delta(x-x_{a})$ for
$z\rightarrow 0$. If we imagine the electromagnetic field propagating in the
medium described by $n(x,z)$ as a bundle of optical rays, then the diffraction
kernel can be interpreted as a path integral over all the possible
trajectories of the optical rays contained in the field as

$\displaystyle K(x_{b},x_{a},z)$ $\displaystyle=$
$\displaystyle\int_{x_{a}}^{x_{b}}\,Dx(z)\,\exp{\left[\frac{i}{\lambdabar}S(x_{b},x_{a},z)\right]}$
(56)

where $\lambdabar=\lambda/2\pi=1/k$, and the action functional $S$ for an
optical ray propagating in the medium described by Eq. (53) is defined as

$S=\int_{0}^{z}L(x,\dot{x},z)dz.$ (57)

The Lagrangian $L$ for an optical ray propagating in a medium with refractive
index $n(x,z)$ can be written as follows

$L(x,\dot{x},z)=n(x,z)\sqrt{1+\dot{x}^{2}}.$ (58)

To justify our starting assumption, _i.e._ , that the electromagnetic field in
the medium can be seen as a collection of rays, we can use different
arguments. One possibility would be to represent the field in its plane wave
components and consider each plane wave as an optical ray Jackson (1986).
Another, more inspiring, possibility is to notice that the paraxial equation
(55) is formally equivalent to the Schrödinger equation for a quantum particle
of mass $n_{0}$ in a potential $n(x,z)$, where the propagation direction $z$
plays the role of time, and $k$ plays the role of $\hbar$. Thanks to this
formal analogy, we can identify a single optical ray as a (massive) photon
propagating in the medium, and then easily understand as the diffraction
kernel can be seen as a path integral over all the possible trajectories that
such quantum particle can take when evolving from the initial state to the
final one.

In the form above, the Lagrangian is of little use, since it provides no
analytical solution for the trajectory of the optical ray and, by extension,
doesn’t really allow for an easy handling of the path integral in Eq. (56). To
circumvent this problem we can however assume that the medium is weakly
inhomogeneous, so that $\Delta n=n(x,z)-n_{0}\ll n$ holds, and we only
consider rays propagating in a small region around the $z$-axis, which
corresponds to assuming $\dot{x}\ll 1$ (this is equivalent to assuming that
the paraxial approximation holds). With these assumptions, we can Taylor
expand the square root in Eq. (58) and the refractive index profile obtaining,
to the leading order in $\dot{x}$ and $x$

$L(x,\dot{x},z)\simeq
n_{0}\left[1+\frac{\dot{x}^{2}}{2}-g^{2}(z)\frac{x^{2}}{2}\right].$ (59)

Notice, that the optical Lagrangian is quadratic in both $x$ and $\dot{x}$,
and can be therefore interpreted as the Lagrangian of a (shifted) harmonic
oscillator with mass $n_{0}$ and $z$-dependent resonance frequency
$\Omega(z)=g(z)$. Thanks to this analogy, we can calculate the kernel
$K(x_{b},x_{a},z)$ in the same manner we did for the harmonic oscillator in
the previous section, with the difference that now we need to account for the
fact that the frequency of the oscillator varies with $z$. In particular, we
can employ the same trick of writing the components of the trajectories as
$x(z)=x_{cl}(z)+y(z)$, with the boundary conditions that the deviation $y(z)$
is zero at the endpoints $z=0$ and $z=z$, and we can then write the
propagator, in analogy with Eq. (39) as

$K(x_{b},x_{a},z)=F(z)\exp{\left[\frac{i}{\lambdabar}S_{cl}(x_{b},x_{a},z)\right]},$
(60)

where

$F(z)=\sqrt{\frac{i}{2\pi\lambdabar}\frac{\partial^{2}}{\partial x_{b}\partial
x_{a}}S_{cl}(x_{a},x_{b},z)}.$ (61)

Following this line of reasoning, and representing the classical action
$S_{cl}(x_{a},x_{b},z)$ in terms of the so-called paraxial ($H_{1}(z)$) and
field ($H_{2}(z)$) rays of a general optical system, we can obtain a similar
result, in our $(1+1)$-dimensional model, than the $(2+1)$-dimensional result
obtained by Gõmez-Reino and J. Liñares Gómez-Reino and Liñares (1987), namely
333Notice that the pre-factor of Eq. (62) differs slightly from Eq. (25) of
Ref. 31. This is due to the different dimensionality of the problem at hand.
While Ref. 31 deals with a $(2+1)$-dimensional problem, in our model we are
only considering one transverse dimension, and therefore the pre-factor has
the form of $\sqrt{X}$ instead of $X$.

$\displaystyle K(x_{b},x_{a},z)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{kn_{0}}{2\pi iH_{1}(z)}}\exp{(ik\,n_{0}\,z)}$ (62)
$\displaystyle\times$
$\displaystyle\exp{\Bigg{[}i\,k\,n_{0}\,\frac{\dot{H}_{1}(z)x_{b}^{2}+H_{2}(z)x_{a}^{2}-2x_{a}x_{b}}{2H_{1}(z)}\Bigg{]}},$

where the dot stands for derivative with respect to $z$ and $H_{1,2}(z)$ are
two independent solutions of the following equation of motion

$\ddot{H}_{1,2}+g^{2}(z)\,H_{1,2}=0.$ (63)

This result, however, is not particularly insightful, and the solution
presented above (or its two-dimensional counterpart presented in Ref. 31) is
quite hard to intuitively link to known results. For this reason, we present
below a much clearer and intuitive approach, which we hope will help the
reader in appreciating the universality and versatility of path integrals
beyond quantum mechanics.

###  IV.2 Paraxial Propagation as a Harmonic Oscillator

Let us go back to the analogy to the Lagrangian (59) and calculate the path
integral deriving from it. The fact that the frequency of the oscillator is
now $z$-dependent doesn’t really allow us to repeat one-to-one the
calculations in Sect. III.4. In particular, after we introduce the deviations
$y(z)$, we cannot represent $F(z)$ in terms of a Fourier series anymore, since
now the oscillator frequency depends on $z$ as well. Instead of doing that,
then, we just follow the same line of reasoning that we used to calculate the
path integral for a free particle in Sect. III.1, namely we divide the
propagation “interval" $Z=z_{b}-z_{a}$ into $N$ smaller intervals of length
$\varepsilon=z_{i+i}-z_{i}$, such that $Z=N\varepsilon$, so that we can write
the action correspondent to the Lagrangian (59) as

$S_{i}(x_{i},x_{i-1})=\frac{n_{0}}{2\varepsilon}(x_{i}-x_{i-1})^{2}-\frac{n_{0}\varepsilon}{2}\Omega_{i}^{2}x_{i}^{2},$
(64)

where $x_{0}\equiv x_{a}$, $x_{N}\equiv x_{b}$, and $\Omega_{i}^{2}\equiv
g^{2}(z_{i})$. This allows us to approximate the path integral in Eq. (56)
using Eq. (24), and evaluate it first over a finite set of $N$ trajectories,
and then to get back to the actual result by taking the limit
$N\rightarrow\infty$. By doing this, moreover, we gain the advantage, that in
each infinitesimal interval $\varepsilon$, the frequency $\Omega_{i}$ of the
oscillator is constant, and so we can use the results for the standard
oscillator within each interval of length $\varepsilon$. If we then introduce
the quantities

$\displaystyle\beta$ $\displaystyle=\frac{n_{0}}{\varepsilon\lambdabar},$
(65a) $\displaystyle\alpha_{i}$
$\displaystyle=\beta\left(1-\frac{1}{2}\Omega_{i}^{2}\varepsilon^{2}\right),$
(65b)

we can then rewrite the path integral in Eq. (56) as $N$ nested Gaussian
integrals, _i.e._ ,

$\displaystyle K_{N}(x_{b},x_{a},z)$ $\displaystyle=$
$\displaystyle\exp{\left[\frac{i\beta}{2}\left(x_{a}^{2}+x_{b}^{2}\right)\right]}$
(66) $\displaystyle\times$ $\displaystyle\int\,dx_{1}\cdots
dx_{N-1}\,\exp{\left[i\sum_{n=1}^{N-1}\alpha_{n}x_{n}^{2}\right]}$
$\displaystyle\times$ $\displaystyle\exp{\left(-i\beta
x_{0}x_{1}\right)}\cdots\exp{\left(-i\beta x_{N-1}x_{N}\right)},$

and then the diffraction kernel can be calculated as

$K(x_{b},x_{a},z)=\lim_{N\rightarrow\infty}K_{N}(x_{b},x_{a},z).$ (67)

The explicit expression for $K_{N}$ can be calculated by recursively applying
the following Gaussian integral result to Eq. (66)

$\displaystyle I$ $\displaystyle=$
$\displaystyle\int\,dx\,\exp{\left\\{i\left[\alpha
x^{2}-\left(a+b\right)x\right]\right\\}}$ (68) $\displaystyle=$
$\displaystyle\sqrt{\frac{i\pi}{\alpha}}\exp{\left[-i\frac{a^{2}+b^{2}}{4\alpha}\right]}\exp{\left[-i\frac{ab}{2\alpha}\right]}.$

This procedure will lead to the quite compact result

$K_{N}(x_{b},x_{a},z)=\sqrt{\frac{a_{N}}{2\pi}}\exp{\left[i\left(p_{N}x_{a}^{2}+q_{N}x_{b}^{2}\right)-a_{N}x_{a}x_{b}\right]},$
(69)

where $a_{N}$, $p_{N}$, and $q_{N}$ are quantities that depend on $\beta$ and
$\alpha_{i}$, and their explicit expression can be found in Ref. 97. We can
then use the result above and take its limit for $N\rightarrow\infty$ to get
the final form of the propagator, which is given explicitly as follows

$\displaystyle K(x_{b},x_{a},z)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{n_{0}\sqrt{\dot{\gamma}_{a}\dot{\gamma}_{b}}}{2\pi
i\lambdabar\sin\phi(z)}}\exp{\left[\frac{in_{0}}{2\lambdabar}\left(\frac{\dot{s}_{b}}{s_{b}}x_{b}^{2}-\frac{\dot{s}_{a}}{s_{a}}x_{a}^{2}\right)\right]}$
(70) $\displaystyle\times$
$\displaystyle\exp{\Bigg{\\{}\frac{in_{0}}{2\lambdabar\sin\phi(z)}\Bigg{[}(\dot{\gamma}_{b}x_{b}^{2}+\dot{\gamma}_{a}x_{a}^{2})\cos\phi(z)}$
$\displaystyle-$ $\displaystyle
2\sqrt{\dot{\gamma}_{b}\dot{\gamma}_{a}}\,x_{a}x_{b}\Bigg{]}\Bigg{\\}}$

where the proper limit for the quantities $a_{N}$, $p_{N}$, and $q_{N}$ has
been taken as instructed in Ref. 97 (also see the same discussion for the
simple harmonic oscillator in Sect. III.4), $\phi(z)=\gamma(z)-\gamma(0)$, and
the quantities $s(z)$ and $\gamma(z)$ are determined from the solution of the
following differential equation

$\frac{d^{2}}{dz^{2}}\xi(z)+\Omega^{2}(z)\xi(z)=0$ (71)

where $\xi(z)=s(z)\exp{[i\gamma(z)]}$, and $\Omega(z)\in\mathbb{R}$. First of
all, notice the similarity between this result and the form of the propagator
for a harmonic oscillator, as given by Eq. (52). In fact, the result above
reduces to the propagator of a harmonic oscillator with constant frequency
$\omega$ by means of the substitution $s(z)=\sqrt{n_{0}/\omega}$ and
$\gamma(z)=\omega\,z$. We can then interpret the propagation of light in a
medium described by the refractive index (53) as basically being given by the
propagator of a harmonic oscillator with a suitably chosen $z$-dependent
frequency (which depends on the longitudinal properties of the refractive
index). The propagator, moreover, is uniquely determined by the amplitude and
phase of the solution of a harmonic-oscillator-like equation of motion for
$\xi(z)$, where the $z$-dependent profile of the medium determines, again, the
characteristic frequency of the oscillator.

Following Ref. 97 we can make this connection appear more evident, by
rewriting the propagator above in terms of the eigenstates of the refractive
index potential given by Eq. (53). This can be done by first rewriting the
$\sin\phi(z)$ and $\cos\phi(z)$ as complex exponentials and then using
Mehler’s formula Olver  _et al._ (2010)

$\displaystyle\frac{\exp{\left[-(X^{2}+Y^{2}-2X\,Y\,Z)/(1-Z^{2})\right]}}{\sqrt{1-Z^{2}}}$
(72) $\displaystyle=$
$\displaystyle\exp{\left[-\left(X^{2}+Y^{2}\right)\right]}\sum_{n=0}^{\infty}\frac{Z^{n}}{2^{n}n!}\text{H}_{n}(X)\text{H}_{n}(Y)$

with $Z=\exp{(-i\phi(z))}$, $X=x_{b}\sqrt{n_{0}\dot{\gamma}_{b}/\lambdabar}$,
and $Y=x_{a}\sqrt{n_{0}\dot{\gamma}_{a}/\lambdabar}$ to transform the
resulting expression in terms of Hermite polynomials and finally obtain (we
generalise $z=0$ to a generic $z=z_{a}$ for convenience)

$K(x_{b},x_{a},z_{a},z_{b})=\sum_{n=0}^{\infty}\psi_{n}^{*}(x_{a},z_{a})\psi_{n}(x_{b},z_{b}),$
(73)

where $\psi_{n}(x,z)$ are the eigenfunctions of a harmonic oscillator with
$z$-dependent frequency $\Omega(z)$, and their explicit expression is given as

$\displaystyle\psi_{n}(x,z)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{1}{2^{n}n!}\sqrt{\frac{n_{0}\dot{\gamma}}{\pi\lambdabar}}}\exp{\left[-i\left(n+\frac{1}{2}\right)\gamma(z)\right]}$
(74) $\displaystyle\times$
$\displaystyle\exp{\left[\frac{in_{0}}{2\lambdabar}\left(\frac{\dot{s}}{s}+i\dot{\gamma}\right)x^{2}\right]}\text{H}_{n}\left(\sqrt{\frac{n_{0}\dot{\gamma}}{\lambdabar}}x\right),$

where $s(z)$ and $\gamma(z)$ are solutions of Eq. (71). For light propagating
in vacuum, $g(z)=0$ (and, therefore $\Omega(z)=0$) and the above expression of
the propagator reduces to the well-known result of the resolution of the
diffraction kernel in terms of Hermite-Gaussian eigenstates of the paraxial
equation Steuernagel (2005). The expression above is moreover fully equivalent
to that found by C. Gõmez-Reino and J. Liñares Gómez-Reino and Liñares (1987),
but it gives more insight on how path integrals can be used to easily solve
complicated problems in optics, such as the propagation of light in a GRIN
medium, whose analytical solution, even in the simplest cases, in not really
available. By establishing the analogy with a time-dependent harmonic
oscillator, on the other hand, path integrals give a rather elegant and
remarkably simple result, which can be intuitively hinted at using simple
arguments, such as the propagation of light in free space.

##  V An Example from Quantum Optics: Path Integral Description of Degenerate
Parametric Amplifiers

Path integrals, in the form defined by Eq. (9), have also been used to
approach various problems in nonlinear optics. In this case, as we will see
throughout this section, the integral over all the possible trajectories of
the quantum particle will be replaced, in the formulation originally proposed
by Hillery and Zubairy in 1982 Hillery and Zubairy (1982), by an integral over
all possible configurations in the complex plane spanned by coherent states.

As an explicit example, we consider here the case of parametric amplification,
whereby a signal incident on an optically nonlinear material is amplified Fox
(2006); Boyd (2008). A pump field impinging on a material characterised by a
$\chi_{2}$ nonlinearity (see Fig. 6) can combine with a signal field, leading
to a depleted pump and an amplified signal (and an idler field, in order to
conserve energy).

In this section, the contents of which based on Ref. 38, we will then apply
the tools presented so far in this tutorial to the case of the degenerate
parametric amplifier, thus providing an example of how path integrals can
prove useful to solve problems in quantum optics.

Figure 6: Pictorial representation of optical parametric amplification. A
nonlinear material with second-order susceptibility ($\chi_{2}$) is impinged
upon by pump photons at angular frequency $\omega_{\rm{p}}$ and signal photons
at $\omega_{\rm{s}}$. Due to the excitation of the medium to a virtual energy
level by the pump (shown at the bottom of the figure), the signal photon
stimulates the emission of photons at $\omega_{\rm{s}}$ and, due to energy
conservation, produces an idler field at $\omega_{\rm{i}}$. For the case of a
degenerate parametric amplifier, the signal and idler fields have the same
frequency, _i.e._ , $\omega_{\rm{s}}=\omega_{\rm{i}}$.

For the sake of simplicity, let us consider a single mode of the radiation
field (generalisations of this formalism to multimode fields is covered
explicitly in Ref. 38), and we represent the field using coherent states, a
natural choice here as the Hamiltonians under consideration will be expressed
in terms of creation and annihilation operators ($a^{\dagger}$ and $a$,
respectively). This assumption, in particular, will allow us to define the
path integral in terms of coherent states, _i.e._ , as a path integral in the
field’s phase space.

###  V.1 Propagator as Path Integral Over Coherent States

We begin by defining the time-evolution operator as $\hat{U}(t_{b},t_{a})$,
which evolves the system’s state at time $t_{a}$ to the state at time $t_{b}$
by $|\psi(t_{b})\rangle=\hat{U}(t_{b},t_{a})|\psi(t_{a})\rangle$ Sakurai
(1994). Denoting the time-dependent Hamiltonian of the system by $H(t)$, the
time-evolution operator is then
$\hat{U}(t_{b},t_{a})=\hat{T}\mathrm{exp}\left[-i\int_{t_{a}}^{t_{b}}H(t^{\prime})dt^{\prime}\right]$,
where $\hat{T}$ is the time-ordering operator Lancaster and Blundell (2014).

If the electromagnetic field is represented in terms of coherent states, we
can readily give a definition of the propagator using the time-evolution
operator defined above as follows

$K(\alpha_{b},t_{b};\alpha_{a},t_{a})=\langle\alpha_{b}|\hat{U}(t_{b},t_{a})|\alpha_{a}\rangle,$
(75)

where $|\alpha\rangle$ represent the set of coherent states, defined as the
eigenstates of the annihilation operator, _i.e._ ,
$\hat{a}|\alpha\rangle=\alpha|\alpha\rangle$ Gerry and Knight (2005). If we
introduce the notation $|\alpha,t\rangle\equiv\hat{U}(t,0)|\alpha\rangle$, we
can rewrite the expression above in the following form, which will prove
useful in the remainder of this section

$K(\alpha_{b},t_{b};\alpha_{a},t_{a})=\langle\alpha_{b},t_{b}|\alpha_{a},t_{a}\rangle.$
(76)

To understand why the form above is useful to our means, let us show how the
quantity above naturally appears when calculating expectation values of
operators and, more generally, correlation functions, in the so-called
$P$-representation. We assume that at time $t=0$ the density matrix of the
field can be written as

$\rho=\int d^{2}\alpha P(\alpha)|\alpha\rangle\langle\alpha|,$ (77)

where $P(\alpha)$ is the Glauber—Sudarshan $P$-function Sudarshan (1963);
Glauber (1963). In this representation, the expectation value of any operator
$\hat{O}(t)$ in the Heisenberg picture can then be written as

$\langle\hat{O}(t)\rangle=\mathrm{Tr}[\rho\hat{O}(t)]=\int d^{2}\alpha
P(\alpha)\langle\alpha|\hat{O}(t)|\alpha\rangle.$ (78)

Using the expression above in combination with the completeness relation of
coherent states, _i.e._ ,

$\frac{1}{\pi}\,\int d^{2}\alpha\,|\alpha,t\rangle\langle\alpha,t|=1,$ (79)

where $d^{2}\alpha=d\operatorname{Re}\alpha\,d\operatorname{Im}\alpha$ spans
the complex-plane defined by coherent states Gerry and Knight (2005), any
normal-ordered correlation function of the form
$\langle\hat{a}^{\dagger}(t_{1})\hat{a}^{\dagger}(t_{2})\cdots\hat{a}(t_{N-1})\hat{a}(t_{N})\rangle$
can be then expressed in terms of the propagator (76). For the simple case of
$\langle\hat{a}(t)\rangle$, this can be readily shown as follows

$\displaystyle\langle\hat{a}(t)\rangle$ $\displaystyle=$
$\displaystyle\int\,d^{2}\alpha\,P(\alpha)\langle\alpha|\hat{a}(t)|\alpha\rangle$
(80) $\displaystyle=$
$\displaystyle\int\,d^{2}\alpha\,P(\alpha)\langle\alpha|\hat{U}^{-1}(t,0)\,\hat{a}\,\hat{U}(t,0)|\alpha\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{\pi}\,\int\,d^{2}\alpha\,d^{2}\,\beta\,P(\alpha)\,\beta\,\left|K(\beta,t;\alpha,0)\right|^{2},$

where to pass from the second to the third line we have employed the
completeness relation for coherent states. For correlation functions
containing more creation and annihilation operators, more propagators,
calculated at different times, will appear in the expression above Hillery and
Zubairy (1982). The propagator is then clearly an important quantity, to
evaluate such expectation values.

We now show how to write the propagator as a path integral. To start with, let
us assume the Hamiltonian of our system is normal-ordered, _i.e._ ,
$H=H(\hat{a}^{\dagger},\hat{a},t)$, and as we did in Sect. III.1 assume to
divide the evolution interval $t_{b}-t_{a}$ into $N$ slices of length
$\varepsilon=(t_{b}-t_{a})/N$. We can mirror this choice directly into Eq.
(76) by inserting $N$ times the completeness relation (79) and set
$t\rightarrow t_{j}=t_{j-1}+\varepsilon$ for each entry. If we do so, we then
obtain

$\displaystyle K(\alpha_{b},t_{b};\alpha_{a},t_{a})$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\pi}\right)^{N}\int d^{2}\alpha_{1}\ldots
d^{2}\,\alpha_{N}\langle\alpha_{b},t_{b}|\alpha_{N},t_{N}\rangle$ (81)
$\displaystyle\times$
$\displaystyle\langle\alpha_{N},t_{N}|\alpha_{N-1},t_{N-1}\rangle\ldots\langle\alpha_{1},t_{1}|\alpha_{a},t_{a}\rangle.$

The quantities $\langle\alpha_{j},t_{j}|\alpha_{j-1},t_{j-1}\rangle$ can be
readily evaluated using the time-evolution operator defined above, and
noticing that since $t_{j}-t_{j-1}=\varepsilon\ll 1$, we can Taylor expand the
time-evolution operator to obtain

$\hat{U}(t_{j},t_{j-1})\simeq
1-i\,\varepsilon\,\int_{t_{j-1}}^{t_{j}}\,d\tau\,H(\hat{a}^{\dagger},\hat{a},\tau).$
(82)

With this result at hand, the individual terms
$\langle\alpha_{j},t_{j}|\alpha_{j-1},t_{j-1}\rangle$ can be brought, after a
simple algebraic manipulation Hillery and Zubairy (1982), in the following
form

$\displaystyle\langle\alpha_{j},t_{j}|\alpha_{j-1},t_{j-1}\rangle$
$\displaystyle\approx$
$\displaystyle\exp{\left[(-1/2)\left(|\alpha_{j}|^{2}+|\alpha_{j-1}|^{2}\right)+\alpha_{j}^{*}\alpha_{j-1}\right]}$
(83) $\displaystyle\times$ $\displaystyle\exp{\left[-i\epsilon
H(\alpha_{j}^{*},\alpha_{j-1},t_{j-1})\right]},$

where we have defined

$H(\beta^{*},\alpha,t)=\frac{\langle\beta|H(a^{\dagger},a,t)|\alpha\rangle}{\langle\beta|\alpha\rangle}.$
(84)

Now that we have an expression of the propagator in terms of $N$ discrete
“trajectories" (_i.e._ , $N$ different coherent states $|\alpha_{j}\rangle$),
we can take the limit $N\rightarrow\infty$ and arrive at the definition of the
path integral in coherent state representation, following the same line of
reasoning used in Sect. III. The details of this calculation are reported in
Ref. 38, and we refer the interested reader therein. The final result of this
calculation is then given as follows:

$\displaystyle K(\alpha_{b},t_{b};\alpha_{a},t_{a})$ $\displaystyle=$
$\displaystyle\int\mathcal{D}\alpha(\tau)\exp{\Bigg{\\{}\int_{t_{a}}^{t_{b}}d\tau\Big{[}\frac{\alpha\dot{\alpha}^{*}-\alpha^{*}\dot{\alpha}}{2}}$
(85) $\displaystyle-$ $\displaystyle
iH(\alpha^{*},\alpha;\tau)\Big{]}\Bigg{\\}},$

with the integration measure taken to mean an integral over all coherent
states parameterised by $\tau$, with $\alpha(t_{a})\equiv\alpha_{a}$ and
$\alpha(t_{b})\equiv\alpha_{b}$.

###  V.2 Propagator for Quadratic Hamiltonians

We now apply the above results, focusing on the class of Hamiltonians which
are at most quadratic in the creation and annihilation operators. This class
of Hamiltonians is of particular interest in quantum optics, since it
describes second-order nonlinear phenomena, within the so-called undepleted
pump approximation Boyd (2008), and it can sometimes also be used to describe
harmonic generation in third-order nonlinear systems Boyd (2008); Drummond and
Hillery (2014).

The most general quadratic Hamiltonian can be written as

$H(\hat{a}^{\dagger},\hat{a},t)=\omega(t)\hat{a}^{\dagger}\hat{a}+\left[f(t)\hat{a}^{2}+g(t)\hat{a}+\text{h.c.}\right],$
(86)

with $f(t)$ and $g(t)$ arbitrary, but well-behaved, time-dependent functions.
The reader might notice a similarity between the Hamiltonian above and that of
a forced harmonic oscillator Feynman and Hibbs (1965). For this class of
Hamiltonians, the integrals appearing in Eq. (85) are all Gaussian, and can be
performed using the methods described in Sect. III. In particular, to
calculate the path integral with the Hamiltonian defined above, one first
needs to discretise the paths over the coherent states $\alpha$ (in the same
manner we discretised the trajectories for the harmonic oscillator in Sect.
III), then express the exponential appearing in Eq. (85) in terms of
$\operatorname{Re}\\{\alpha\\}$ and $\operatorname{Im}\\{\alpha\\}$
explicitly, and notice, that the correspondent integrals are Gaussian in both
the real and the imaginary parts of $\alpha$. After having calculated a single
term, one could then calculate the remaining integrals iteratively, as we have
done for the examples in Sect. III. After a lengthy but straightforward
calculation, which is partially covered in Appendix A of Ref. 38, the
propagator for a quadratic Hamiltonian assumes the following form

$K(\bm{\alpha},\textbf{t})=F(\bm{\alpha},\textbf{t})\exp{\left[-i\,\Sigma(\bm{\alpha},\textbf{t})\right]},$
(87)

where $\bm{\alpha}=\\{\alpha_{a},\alpha_{b}\\}$,
$\textbf{t}=\\{t_{a},t_{b}\\}$ and the functions $F$ and $\Sigma$ are defined
as

$\displaystyle F(\bm{\alpha},\textbf{t})$ $\displaystyle=$
$\displaystyle\exp{\Big{[}-\frac{1}{2}(|\alpha_{b}|^{2}+|\alpha_{a}|^{2})+Y(t_{b})\alpha_{b}^{*}\alpha_{a}}$
(88) $\displaystyle+$ $\displaystyle
X(t_{b})(\alpha_{b}^{*})^{2}+Z(t_{b})\alpha_{b}^{*}\Big{]},$
$\displaystyle\Sigma(\bm{\alpha},\textbf{t})$ $\displaystyle=$
$\displaystyle\int_{t_{a}}^{t_{b}}\,d\tau\Big{\\{}f(\tau)\Big{[}2X(\tau)+Z^{2}(\tau)+\alpha_{a}^{2}Y^{2}(\tau)$
(89) $\displaystyle+$ $\displaystyle
2\alpha_{a}Y(\tau)Z(\tau)\Big{]}+g(\tau)\Big{[}Z(\tau)+\alpha_{a}\,Y(\tau)\Big{]}\Big{\\}},$

where the auxiliary function $X(t)$ is constrained by

$\frac{dX}{dt}=-2i\omega(t)X-4if(t)X^{2}-if^{*}(t),$ (90)

with initial condition $X(t_{a})=0$. The functions $Y(t)$ and $Z(t)$ are
instead defined in terms of $X(t)$ as

$Y(t)=\mathrm{exp}\left(-i\int_{t_{a}}^{t}d\tau\left[\omega(\tau)+4f(\tau)X(\tau)\right]\right),$
(91)

and

$\displaystyle Z(t)$ $\displaystyle=$
$\displaystyle-i\int_{t_{a}}^{t}d\tau\left[g^{*}(\tau)+2g^{*}(\tau)X(\tau)\right]$
(92) $\displaystyle\times$
$\displaystyle\exp{\left(-i\int_{\tau}^{t}d\tau^{\prime}\left[\omega(\tau^{\prime})+4f(\tau^{\prime})X(\tau^{\prime})\right]\right)}.$

###  V.3 Propagator for Degenerate Parametric Amplification

A degenerate parametric amplifier is characterised by the following quadratic
Hamiltonian Boyd (2008)

$H(t)=\omega a^{\dagger}a+\kappa\left(e^{2i\omega t}a^{2}+e^{-2i\omega
t}a^{\dagger 2}\right),$ (93)

where $\omega$ is the angular frequency of the mode and $\kappa$ is a coupling
constant. The above Hamiltonian can clearly be seen to fall into the class of
Hamiltonians given in Eq. (86) if we identify $\omega(t)=\omega$, $f(t)=\kappa
e^{2i\omega t}$, and $g(t)=0$. Although parametric amplification is formally a
3-wave process involving, as depicted in Fig. 6, a pump, a signal, and an
idler field, it is common practice in nonlinear optics experiments to work
within the so-called undepleted pump approximation Boyd (2008); Drummond and
Hillery (2014), which treats the pump mode as a classical (bright) field,
whose number of photons does not change significantly (_i.e._ , the pump field
remains undepleted) during the nonlinear interaction. Within this
approximation, then, the Hamiltonian becomes effectively quadratic in the
signal and idler modes, and can be described by Eq. (93) above. The explicit
form of the propagator given below is then valid within this approximation. A
fully quantised analysis of parametric amplification beyond the undepleted
pump approximation in terms of path integrals has been presented by Hillery
and Zubairy in Ref. 37.

In this case, with the form of $\omega(t)$, $f(t)$, and $g(t)$ given above,
the functions $X(t)$, $Y(t)$, and $Z(t)$ can be written in an explicit form
(in particular, Eq. (90) cabn be solved analytically), giving $Z(t)=0$ and
Hillery and Zubairy (1982)

$\displaystyle X(t)$ $\displaystyle=\frac{1}{2i}e^{-2i\omega
t}\mathrm{tanh}\left[2\kappa(t-t_{a})\right],$ (94a) $\displaystyle Y(t)$
$\displaystyle=e^{-i\omega(t-t_{i})}\mathrm{sech}\left[2\kappa(t-t_{a})\right].$
(94b)

These time-dependent quantities can then be directly substituted into Eq.
(87)), finally yielding the propagator for the degenerate parametric
amplifier:

$\displaystyle K(\alpha_{b},t_{b};\alpha_{a},t_{a})$ $\displaystyle=$
$\displaystyle\sqrt{\mathrm{sech}\left[2\kappa(t_{b}-t_{a})\right]}\exp{\left(-\frac{|\alpha_{b}|^{2}+|\alpha_{a}|^{2}}{2}\right)}$
(95) $\displaystyle\times$
$\displaystyle\exp{\Bigg{\\{}\alpha_{b}^{*}\alpha_{a}e^{-i\omega(t_{b}-t_{a})}\mathrm{sech}\left[2\kappa(t_{b}-t_{a})\right]}$
$\displaystyle-$ $\displaystyle
i\left(\frac{\alpha_{b}^{*}}{\sqrt{2}}\right)^{2}e^{-2i\omega
t_{b}}\mathrm{tanh}\left[2\kappa(t_{b}-t_{a})\right]$ $\displaystyle-$
$\displaystyle i\left(\frac{\alpha_{a}}{\sqrt{2}}\right)^{2}e^{2i\omega
t_{a}}\mathrm{tanh}\left[2\kappa(t_{b}-t_{a})\right]\Big{]}\Bigg{\\}}.$

##  VI Predicting Optical Properties of Matter Using Path Integrals

The examples discussed above only deal with a single particle, either freely
propagating, or interacting with a suitable potential. In all these cases, the
complexity of the problem is low enough, to allow analytical solutions to the
path integrals. As soon as many-body effects are taken into account, as it is
the case, for example, for atoms and molecules, handling the path integral
analytically becomes impossible, and numerical techniques for its evaluation
are necessary. A particularly successful computational method for this task
has been Path Integral Monte Carlo (PIMC). The interested reader can download
David Ceperley’s group PIMC++ open-source code for Path Integral Monte Carlo
simulations, available, together with its relative documentation, at
https://pimc.soft112.com.

In this section, we then present some of the results obtained using PIMC to
predict the _exact_ linear and nonlinear optical properties of simple systems
as a function of frequency. These results, recently demonstrated by Tiihonen
_et al._ Tiihonen (2019), constitute one of the first attempts to use PIMC to
investigate the optical properties of matter. There, however, only very simple
systems were studied, such as the hydrogen (H) and hydrogen-like
($\text{Li}^{+}$ and $\text{Be}_{2}^{+}$) atoms, the helium atoms He and
$\text{He}^{+}$, the hydrogen molecule $\text{H}_{2}$ and $\text{H}_{2}^{+}$,
hydrogen-helium ($\text{HeH}^{+}$) and hydrogen-deuterium ($\text{HD}^{+}$)
molecules, and the positronium atom. These systems have been studied employing
different methods, including finite-field simulations for static
polarisabilities Tiihonen, Kylänpää, and Rantala (2015), polarisability
estimators for simulation without the external field Tiihonen, Kylänpää, and
Rantala (2016), static field-gradient polarisabilities Tiihonen, Kylänpää, and
Rantala (2017), and finally, dynamic polarisabilities and van der Waals
coefficients Tiihonen, Kylänpää, and Rantala (2018). The last one, in
particular, is of great interest, as the macroscopic electric susceptibility
is constructed starting from the dynamic frequency dependent polarisabilities.

The results presented in the aforementioned references, and summarised in this
section represent an exceptional benchmark to gauge the capabilities of PIMC
in providing accurate and multidimensional information about the optical
properties of matter, and hint at the potential impact PIMC could have in
photonics, as a new modelling platform to calculate _exactly_ the optical
response of exotic materials, such as 2D materials, the understanding of which
is still in its infancy Anton  _et al._ (2018). A great advantage PIMC would
give, for example, is the possibility to investigate the nonlinear properties
of 2D materials far from equilibrium, a physical regime that is currently
poorly understood. This could be of particular importance for the nonlinear
optical properties of 2D materials, since for such materials, non-equilibrium
dynamics can be easily reached at relatively small optical powers.

###  VI.1 Polarisability in PIMC

The optical response of matter Jackson (1986), frequently described in terms
of the matter polarisation vector P, is completely determined by electron
dynamics solely within the limits of the the Born-Oppenheimer (BO)
approximation, which specifies that the nuclear dynamics occurs on a much
longer timescale and can be therefore neglected, it is implicitly assumed
Fetter and Walecka (2003). In general, however, when thermal effects are
explicitly taken into account, this approximation might not be valid anymore,
as the role of the nuclei becomes more prominent as the temperature of the
system is increased. To properly take into account these effects, then, an
approach able to go beyond BO is needed to get the correct results. PIMC,
then, is the most viable, if not the only possible, approach to systematically
study the optical properties of materials in this particular regime (which,
for example, includes the calculation of such properties at room temperature).

While optics makes extensive use of the (linear and nonlinear) susceptibility
tensor to describe the optical properties of a material Boyd (2008), from a
PIMC perspective, it is better to work with polarisabilities, as they can be
defined quite easily in terms of path integrals. For example, the
$ij$-component of the linear (_i.e._ , dipole) polarisability tensor, namely
$\alpha_{ij}$, can be simply calculated as the (functional) derivative of the
expectation value (_i.e._ , first order correlator) of the dipole moment
$\bm{\mu}$ oriented, say, along the $i$-direction, with respect to an external
electric field oriented along the $j$ direction Tiihonen (2019), _i.e._ ,

$\alpha_{ij}=\frac{\partial}{\partial E_{j}}\langle\mu_{i}\rangle,$ (96)

where the expectation value of an operator $\hat{O}$ is defined with respect
to the density matrix operator of the system in the so-called imaginary time
representation Tiihonen (2019), _i.e._ ,

$\langle\hat{O}\rangle=\mathcal{N}\operatorname{Tr}\left\\{\hat{O}\,\exp{\left[-\frac{S(T)}{\hbar}\right]}\right\\},$
(97)

where $\mathcal{N}$ is a suitable normalisation constant,
$S(T)=\left(\hat{H}_{0}-\hat{\textbf{Q}}\cdot\textbf{E}\right)/k_{B}T$ is the
temperature-dependent action of the system (described by the Hamiltonian
$\hat{H}_{0}$) interacting with the electromagnetic field E, and
$\hat{\textbf{Q}}$ is the electric moment operator Tiihonen (2019). Notice
that to represent the expectation value above as a path integral, it is
computationally more convenient to transform the complex phase factor
$\exp{\left(i\,S/\hbar\right)}$, naturally arising from path integrals, into
an exponentially decaying term $\exp{\left(-S/\hbar\right)}$ by means of a
Wick rotation Srednicki (2007) ( _i.e._ , a change in reference frame from
real to imaginary time, namely $t\rightarrow i\tau$), since computationally it
is much easier to deal with real-valued, exponentially decaying terms, rather
than complex-valued, spuriously oscillating ones. In fact, the everywhere
positive exponential function can be considered as a probability distribution
for sampling the imaginary time paths with Metropolis Monte Carlo algorithm in
canonical (NVT) ensemble.

PIMC makes use of several different strategies to compute the quantities
defined above. The zero-field polarisability estimators, for example, are
Hellmann—Feynman type operators Tiihonen (2019), whose expectation values give
the polarisabilities of various types and order, without the need of including
an external electric field in the simulation.

Computation of dynamic multipolar polarisabilities, on the other hand, is more
sophisticated and it requires to first determine, through PIMC, the multipole-
multipole correlation function in the so-called imaginary time representation,
then analytically continue it to obtain the spectral function in the real time
domain (a task that is highly non-trivial, due to the fact that it is an ill-
posed numerical problem), and finally both the real and imaginary parts of the
polarisability.

###  VI.2 From Polarisability to Susceptibility

Figure 7: Real (top) and imaginary (bottom) part of the dynamic polarisability
as a function of frequency for Hydrogen (H, left) and Helium (He, right) atoms
at $T=2000$ K, within the BO approximation. In these plots, a blue line
represents the dipole polarisability $\alpha_{1}(\omega)$, while a green line
represents the quadrupole polarisability $\alpha_{2}(\omega)$. The black dots
represent the reference simulations of polarisabilities calculated at $T=0$ K.
Dotted lines correspond to big time steps, i.e., $\Delta\tau=0.1$ for
Hydrogen, and $\Delta\tau=0.025$ for Helium, while solid lines correspond to
small time steps, namely $\Delta\tau=0.05$ for Hydrogen, and
$\Delta\tau=0.0125$ for Helium. In PIMC, time steps are defined, in imaginary
time, as $\Delta\tau=1/(Mk_{B}T)$, where $k_{B}$ is the Boltzmann’s constant,
$T$ is the temperature of the system, and $M$ is the so-called Trotter number,
which controls the accuracy of the discretisation of the propagator
$\exp{(-\Delta\tau\,\hat{H})}$. In these plots, the frequency is measured in
Hartree (1 Ha $=$ 2 Ry $=27.2$ eV). To get a better grasp of the span of the
frequency axis of the plots above, a frequency of $\omega=0.1$ Ha corresponds
to a wavelength $\lambda=456$ nm, while a frequency of $\omega=1$ Ha
corresponds to a wavelength $\lambda=45$ nm. This figure has been taken from
Ref. 107 with permission of the author.

Assuming a very low density gas phase of these small atoms or molecules, one
can evaluate the corresponding macroscopic susceptibility from the atomic and
molecular polarisabilities. At higher densities, interactions and chemical
reactions of these moieties will change the polarisabilities, composition or
both, which consequently leads to the density dependent susceptibilities.

With increasing density and, in particular, in case of liquids and solids,
where several electrons interact strongly, the Fermion Sign Problem (FSP)
emerges. This is the notorious challenge for Monte Carlo evaluation of
essential expectation values of identical fermions, and it discloses partly a
still open problem. The problem emerges from the evaluation of the difference
of two large contributions with opposite signs, resulting from the sign of the
density matrix (or wave function). There are practical solutions to FSP in
PIMC simulations for many cases Ceperley (1996); Pierleoni  _et al._ (1994);
Militzer and Ceperley (2001); Gorelov  _et al._ (2020); Tiihonen (2019), but
not any general robust approaches, yet. Most of the solutions are based on
finding approximate or iterative nodal surfaces of the density matrix.

The real-time path integral (RTPI) approach Tiihonen, Kylänpää, and Rantala
(2016); Ruokosenmäki and Rantala (2015); Ruokosenmäki  _et al._ (2017);
Gholizadehkalkhoran, Ruokosenmäki, and Rantala (2018) may offer a remedy in
the future, as it works directly with the wave function with explicit sign,
even though only at zero Kelvin. However, as for most of the cases of
interests, electrons at room temperature can be well-described with the zero
Kelvin temperature model, this is not an issue. Moreover, including
temperature in RTPI should be, in principle, possible.

Overall, the FSP represents a challenging problem, and sits at the forefront
of current research in PIMC methods. The interested reader can find in Refs.
20; 21; 22; 23, and references therein, a good starting point to delve deeper
in this fascinating, yet challenging, problem. In the examples that we are
going to present below, however, we have purposely selected simple systems,
where there are no more than two electrons in each moiety, where at low
temperatures we can assume the system to be in its singlet ground state, thus
enabling us to directly label the light fermions with opposite spins and
avoiding FSP entirely.

###  VI.3 Some Examples

In this section we present some results concerning simple systems, such as
one- and two-electron atoms and molecules, obtained using PIMC. In particular,
we focus our attention on the electric polarisability, to show what are the
capabilities of this computational method, and how it could be used to provide
a new interesting platform for photonics, to calculate the optical properties
of interesting materials, beyond the standard approaches and approximations
frequently used in optics Jackson (1986); Boyd (2008).

Figure 8: Real (top) and imaginary (bottom) part of the dynamic polarisability
as a function of frequency for molecular Hydrogen $H_{2}$ for various
temperatures. The left panel presents the results for the dipole
polarisability $\alpha_{1}(\omega)$, while the right panel contains the
results for the quadrupole polarisability $\alpha_{2}(\omega)$. The black dots
represent the reference simulations of polarisabilities calculated at $T=0$ K.
Dotted lines correspond to big time steps, while continuous lines correspond
to smaller time steps. For these simulations, the time step has been adapted
according to the chosen temperature, and it varies from $\Delta\tau=0.02$ to
$\Delta\tau=0.05$. In these plots, the frequency is measured in Hartree (1 Ha
$=$ 2 Ry $=27.2$ eV). To get a better grasp of the span of the frequency axis
of the plots above, a frequency of $\omega=0.1$ Ha corresponds to a wavelength
$\lambda=456$ nm, while a frequency of $\omega=1$ Ha corresponds to a
wavelength $\lambda=45$ nm. This figure has been taken from Ref. 107 with
permission of the author.

In both examples below, the polarisability is calculated from the Kubo-like
formula for the time-dependent optical susceptibility Fetter and Walecka
(2003)

$\chi(t)=\frac{i}{\hbar}\Theta(t)\,\langle\,[\hat{P}(t),\hat{Q}(0)]\,\rangle=-G_{ret}(t),$
(98)

where the square brackets denote a commutator, the angled brackets denote
instead thermal average [see Eq. (97)], $\Theta(t)$ is the Heaviside step
function Byron and Fuller (1992), $\hat{Q}$ is a multipole operator (which
represents dipole or quadrupole interactions, for the cases presented below),
$G_{ret}(t)$ is the retarder Green’s function of the perturbation $\hat{P}$
and the latter represents the perturbation generating the optical
susceptibility. The link between susceptibility and polarisability is better
understood in Fourier space, where the following relation is valid Tiihonen,
Kylänpää, and Rantala (2018)

$\chi(\omega)=-\int\,\frac{d\,\Omega}{\pi}\,\frac{A(\Omega)}{\omega-\Omega+i\eta}=\langle\,\alpha(\omega)\,\rangle,$
(99)

where the spectral function $A(\omega)$ is related to the susceptibility via
the Kramers-Kronig relation Jackson (1986), i.e.,
$A(\omega)=\operatorname{Im}\\{\chi(\omega)\\}$. PIMC is then used to
calculate the value of the expectation value in Eq. (98) (or, equivalently,
the retarded Green’s function), and from it both the real and imaginary part
of the polarisability are then calculated. Notice, moreover, that the real
part of the polarisability $\operatorname{Re}\\{\alpha(\omega)\\}$ corresponds
to the actual optical response of the material, while its imaginary part
$\operatorname{Im}\\{\alpha(\omega)\\}$ is, basically the spectral function
$A(\omega)$.

The first example we present is given in Fig. 7, which depicts the real (top)
and imaginary (bottom) part of the dipole (blue line) and quadrupole (green
line) polarisabilities $\alpha_{1}(\omega)$ and $\alpha_{2}(\omega)$ for the
case of atomic Hydrogen (left) and Helium (right) at $T=2000$ K. The
simulations have been performed using PIMC with two distinct characteristic
time steps $\Delta\tau$, one large (solid green and blue lines in Fig. 7) and
one small (dotted green and blue lines in Fig. 7), to rule out possible
numerical artefacts Tiihonen, Kylänpää, and Rantala (2018). The results have
been then compared with polarisability data available in literature, and
obtained using a $T=0$ K approach, with methods other than PIMC Bishop and
Pipin (1993). This reference is indicated by black dots in Fig. 7. For low
frequencies, i.e., in the region $0<\omega<0.4$ Ha (corresponding to
$0<\lambda<1.14$ $\mu m$), PIMC simulations reproduce very well the standard
$T=0$ K result. This, essentially, means, that for the hydrogen atom, the
eigenstates at $T=2000$ K are essentially the same as those calculated at
$T=0$ K, and, as a consequence of that, the low-frequency polarisability in
Fig. 7 fits perfectly the results obtained at zero temperature Bishop and
Pipin (1993).

As the frequency is increased beyond $\omega=0.5$ Ha (corresponding, for the
case of Hydrogen, to its ionisation energy), the zero temperature approach
fails to faithfully describe the optical properties of the system. In this
high frequency range, i.e., ranging from the near infrared to the deep UV,
PIMC produces instead a more accurate and reliable result. This is an
important result. Usually, the models used in photonics to describe the
optical properties of various materials, are based on a zero temperature
approach. Although this is perfectly fine (and enough) for most materials,
since the wavelength at which non-equilibrium effects are starting to become
non-negligible is either in the far IR or the deep UV, this is not the case
for 2D materials, which can be easily be driven out of equilibrium in the
wavelength range of interest for different applications (typically ranging
from THz to visible). In this case, traditional, zero-temperature-based,
methods would fail to correctly describe the influence of out-of-equilibrium
effects onto the optical properties of 2D materials, while PIMC, on the other
hand, are able to exactly account for these contributions, as they naturally
take these into account.

The second example deals with a slightly more complicated system, namely
molecular hydrogen $\text{H}_{2}$. Figure 8 depicts the real (top) and
imaginary (bottom) part of the dipole (left) and quadrupole (right)
polarisabilities $\alpha_{1}(\omega)$ and $\alpha_{2}(\omega)$, respectively,
for several different values of the ambient temperature $T$. A comparison with
reference simulations obtained at $T=0$ K from Ref. 110 is also shown (black
dots in Fig. 8).

Notice, in particular, the deviation from the zero temperature behaviour in
the small frequency regime (insets in the top part of Fig. 8). This
discrepancy is due to the fact, that in the case of molecules, the nuclear
contribution (mainly nuclear vibrations induced by finite temperature) cannot
be neglected, as it significantly changes the behaviour of the polarisability
at low frequencies. This change amounts only to a shift of the value of
$\operatorname{Re}\\{\alpha_{1}(\omega)\\}$ for the dipole polarisability (top
let inset in Fig. 8), while it radically changes the form of
$\operatorname{Re}\\{\alpha_{2}(\omega)\\}$ at low frequencies, for the case
of quadrupole polarisability (top right inset in Fig. 8). Notice, moreover, as
in this case temperature has also a significant effect in shaping the
imaginary part of the quadrupole polarisability at very low frequencies (lower
right inset in Fig. 8).

Although they show the optical properties of very simple systems, such as H,
He and $\text{H}_{2}$, the examples presented in this section offer good
insight on the capabilities and potential of using PIMC for calculating the
optical properties of even more complicated systems. The large frequency range
PIMC can span (from THz to the deep UV) represents an enormous advantage to
test the current models for susceptibilities and polarisabilities of photonic
materials, and could in particular provide tremendous insight on the nonlinear
properties of exotic materials, such as epsilon-near-zero materials, for
example, as it provides a single, _ab-initio_ method to calculate the optical
response in an exact manner, without inserting any approximation or without
introducing different models for different frequency regimes.

Besides that, the examples provided in this section also help the reader in
getting a better grasp of the breadth of path integrals, and how they do not
only represent an elegant, alternative formulation of quantum mechanics, but
also how their use in atomic and molecular physics can be pivotal for the
accurate description of the properties of matter.

##  VII Path Integral for Classical and Quantum Fields

In the previous sections, we have seen how it is possible to describe the
dynamics of a quantum particle (or a collection thereof) in terms of path
integrals. In this section, instead, we extend this formalism to the case of
quantum fields, where the building blocks for the path integral are not
particles and all the trajectories they can possibly take to evolve from the
initial to the final state, but rather all the possible configurations of a
field, that serve the same purpose.

A quantum field can be seen as a continuous collection of harmonic
oscillators, one for each point in space where the field is defined. As such,
then the results presented for the harmonic oscillator in Sect. III.4 can
serve as guidance to better understand this section, and, more generally, the
properties of a QFT. Some of the language and tools used here, however, are
slightly different than those introduced above, and it is worthwhile to
introduce some key concepts and methods proper of QFT before specialising our
attention to the electromagnetic field. To do that in the most beginner-
friendly manner, we use the case of a scalar field as a guidance to construct
the path integral, and show how interactions can be introduced in such a
theory, and what is the consequence of doing that. Then, we will qualitatively
and intuitively extend these results to the case of a vector field, and use
those results to look at the linear and nonlinear dynamics of the
electromagnetic field in an arbitrary medium.

This section follows closely the notation and line of reasoning of Ref. 109.
For a more formal discussion of path integrals in field theory, we invite the
interested reader to consult Refs. 109; 88; 8; 9. For the sake of simplicity,
natural units $c=\hbar=m_{e}=1$ are implicitly assumed. If needed, these
constants can be easily reinserted in the final results of calculations, by
means of simple dimensional analysis. Moreover, summation over repeated
indices (Einstein summation convention) is also understood.

###  VII.1 Path Integral for a Scalar Field

To start with, let us then consider a free scalar field $\phi(x)$ (where $x$
is a shorthand for $\\{\textbf{r},t\\}$), whose Lagrangian density is given by

$\mathcal{L}_{0}(\phi(x))=\frac{1}{2}\eta^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi(x)-\frac{1}{2}m^{2}\phi(x)^{2},$
(100)

where $\eta^{\mu\nu}=\text{diag}(-1,1,1,1)$ is the Minkowski flat space metric
tensor Srednicki (2007), and the subscript $0$ indicated that the Lagrangian
is that of the free field, _i.e._ , $\phi(x)$ is not interacting with either
itself, or its environment.

Before going any further, it is instructive to discuss the physical meaning of
the parameter $m$ appearing in the above equation. In QFT, it is customary to
give $m$ the meaning of mass of the scalar field $\phi$, as it could represent
a scalar massive boson such as the Higgs boson Kaku (1993), or even a scalar
fermion in more complicated theories involving dark matter Kawamura (2014). In
optics, on the other hand, we frequently refer to the case $m=0$, as the
electromagnetic field is composed by photons, which are massless bosons Kaku
(1993). However, when propagating inside a medium, the electromagnetic field
acquires a mass, proportional to the refractive index of the medium Luneburg
(1965). This can be easily proven by recalling that the dispersion relation
for a monochromatic electromagnetic field propagating inside a homogeneous
medium characterised by a refractive index $n$ is given by Jackson (1986)
$k^{2}=\omega^{2}n^{2}$. The dispersion relation for a massive scalar field
$\phi$, on the other hand, can be readily calculated from its equation of
motion (i.e., the Klein-Gordon equation Srednicki (2007)), which can be
readily derived from Eq. (100) using the Euler-Lagrange equations for fields,
i.e.,

$\partial_{\mu}\left[\frac{\partial\mathcal{L}_{0}}{\partial\left(\partial_{\mu}\phi\right)}\right]-\frac{\partial\mathcal{L}_{0}}{\partial\phi},$
(101)

which leads to $(\nabla^{2}-\partial_{t}^{2}-m^{2})\phi=0$. If we now assume
plane wave solutions for $\phi$, the dispersion relation for a massive scalar
field can be written as $k^{2}=\omega^{2}-m^{2}$. Notice, now, that upon
defining $m^{2}=\omega^{2}(1-n^{2})$, the dispersion relation for a massive
scalar field $\phi$ matches that of an electromagnetic field inside a
homogeneous medium. This allows us to identify the mass (squared) of the
scalar field $\phi$ with $m^{2}=1-n^{2}=1-\varepsilon=\chi$, i.e., the
medium’s susceptibility Jackson (1986). In this sense, therefore, while a
massless scalar field describes a (scalar 444A scalar electromagnetic field is
defined as an electromagnetic wave whose polarisation is the same in any point
of its wavefront, and always remains the same during propagation. If the
interaction of a scalar electromagnetic field with a system causes its state
of polarisation to vary, this variation will happen in the same manner in
every point of its wavefront. A good example of a scalar electromagnetic field
is an optical beam, which is typically written as
$\textbf{E}(\textbf{r},t)=\psi(\textbf{r},t)\,\hat{\textbf{u}}$, where
$\hat{\textbf{u}}$ is a constant (unitary) transverse vector (with respect to
the rpopagation direction of the field), which describes the field
polarisation. The dynamics of an optical beams are then fully determined by
$\psi(\textbf{r},t)$, which, in general, is a solution of the scalar wave
equation.) electromagnetic field in vacuum, a massive one describes instead
the (scalar) electromagnetic field in a homogeneous medium.

We now proceed to the definition of path integral for the scalar field $\phi$.
To do so, let us notice that the Lagrangian density in Eq. (100) has a similar
structure to that of a harmonic oscillator, as given in Eq. (36). The first
term in Eq. (100), in fact, is quadratic in the derivative of the field, and
can be associated with its “kinetic energy", while the second term, quadratic
in the field, is its “potential energy". Analogously to what we have done in
the previous sections, we can introduce the action of the field as

$S_{0}(\phi(x),J)=\int\,d^{4}x\,\left[\mathcal{L}_{0}(\phi(x))+J(x)\phi(x)\right],$
(102)

where $d^{4}x=dt\,d^{3}r$ and we have introduced the source term $J(x)$
(analogous to the external force acting on a harmonic oscillator Feynman and
Hibbs (1965)), which, for the sake of this tutorial, can be imagined to simply
be a convenient mathematical tool, introduced to make calculations easier. The
physical meaning of $J(x)$, as its name suggests, is to represent a source of
the field. The $x$-dependence of the field $\phi(x)$ is from here henceforth,
and if not specified otherwise, implicitly assumed, i.e., we define
$\phi\equiv\phi(x)$, to make our notation less cumbersome.

We can then define the path integral for a free scalar field in the same
manner as we did in Eq. (9), thus obtaining

$Z_{0}(J)=\int\,\mathcal{D}\phi\,e^{iS_{0}[\phi,J]}$ (103)

where the functional measure $\mathcal{D}\phi(x)$ plays the same role as the
sum over all possible trajectories $Dx(t)$ in Eq. (9), with the conceptual
difference, however, that in this case the sum (integral) is extended to all
possible configurations of the function $\phi(x)$. Practically, this means
that instead of slicing time into small intervals (as we did in the previous
sections for the case of quantum particles), we divide the volume of spacetime
$\mathcal{V}$ into $N$ elementary cubes, i.e.,
$\mathcal{V}=\bigcup_{k=1}^{N}V_{k}$, with $V_{k}=\delta^{4}$, and assume that
in each $V_{k}$ the field $\phi$ is constant, i.e.,
$\phi\simeq\phi(t_{i},x_{j},y_{n},z_{m})=\phi(x_{k})\equiv\phi_{k}$ (this
implies that since each field $\phi$ is a function of 4 variables, each taking
$N$ possible values, this discretisation of spacetime results in having a
total of $N^{4}$ fields). This then means, that the integration measure can be
written as

$\mathcal{D}\phi=\lim_{N\rightarrow\infty}\prod_{i,j,n,m=1}^{N}d\phi(t_{i},x_{j},y_{n},z_{m})\equiv\lim_{N\rightarrow\infty}\prod_{k=1}^{N^{4}}d\phi_{k}.$
(104)

In doing so, each field configuration $\phi_{k}$ becomes an independent
integration variable, and the integral $\int\,\mathcal{D}\phi$ reduces to a
usual Riemann/Lebesgue integral in $N^{4}$-dimensions.

The above definition of path integral is valid both for a classical, as well
as a quantum field. For the latter, the field $\phi(x)$ needs to be equipped
with a suitable set of commutation relations, that define the quantum nature
and structure of the algebra behind it. To do that, we can revert back to the
analogy with the harmonic oscillator, and use it to find a straightforward
definition of such quantities. First, we discretise the field $\phi(x)$ and
associate a harmonic oscillator to each point in space. Then, for each
oscillator, we introduce the usual commutation relations between position and
momentum, _i.e._ , $[q_{a}(t),p_{b}(t)]=i\delta_{ab}$, where $\delta_{ab}$
simply indicates that the oscillators are mutually decoupled. If we then take
the continuum limit (to revert back to the field $\phi(x)$) the commutation
relations become those of the field $\phi$ itself with its canonically
conjugated momentum $\Pi=\partial\mathcal{L}/\partial(\partial_{t}\phi)$,
namely

$\displaystyle\left[\phi(\textbf{r},t),\Pi(\textbf{r}^{\prime},t)\right]$
$\displaystyle=i\,\delta(\textbf{r}-\textbf{r}^{\prime}),$ (105a)
$\displaystyle\left[\phi(\textbf{r},t),\phi(\textbf{r}^{\prime},t)\right]$
$\displaystyle=\left[\Pi(\textbf{r},t),\Pi(\textbf{r}^{\prime},t)\right]=0.$
(105b)

Equation (103), moreover, has implicitly been defined such that the time
interval used for the time integration in the action spans the whole real
axis, _i.e._ , $-\infty<t<\infty$. Moreover, we require that $Z_{0}(J)$ is
normalised in such a way that $Z_{0}(0)=1$.

For a classical field, $Z_{0}(J)$ represents its partition function
$Z_{0}(J)$, from which, in analogy with statistical mechanics, it is possible
to derive all the relevant quantities of the classical field, such as energy,
momentum, entropy, and Green’s function.

For a quantum field, it represents the vacuum-to-vacuum correlation function,
which allows to calculate transition probabilities and correlation functions
of field operators with respect to the ground state (or vacuum), which is
typically a stable and well-defined state of the QFT at-hand. In other terms,
$Z_{0}(J)$ for a quantum field only describes a QFT at equilibrium (_i.e._ ,
at zero temperature). To deal with non-equilibrium systems, a further step in
the definition of the path integral is needed, which is beyond the scope of
this tutorial. The interested reader is then referred to Ref. 113 for further
information about non-equilibrium QFTs.

###  VII.2 The Propagator for a Free Scalar Field

To get familiar with the formalism of QFT, in this section we calculate the
propagator for the free scalar field, which constitutes, to a certain extent,
the QFT analogue of the free quantum particle described in Sect. III.1. To
this aim, let us first introduce the Fourier transform of the field as

$\displaystyle\varphi(k)$ $\displaystyle=\int\,d^{4}x\,\phi(x)\,e^{-ikx},$
(106a) $\displaystyle\phi(x)$
$\displaystyle=\int\,\frac{d^{4}k}{(2\pi)^{4}}\,\varphi(k)\,e^{ikx},$ (106b)

where $kx=\textbf{k}\cdot\textbf{r}-\omega t$. If we substitute the expression
of $\phi$ above into Eq. (102) we arrive, after some straightforward
calculations, at the following expression for the action in Fourier space

$\displaystyle S_{0}$ $\displaystyle=$
$\displaystyle\int\,\frac{d^{4}k}{(2\pi)^{4}}\,\Big{[}-\frac{1}{2}\varphi(k)\left(k^{2}+m^{2}\right)\varphi(-k)$
(107) $\displaystyle+$
$\displaystyle\tilde{J}(k)\varphi(-k)+\tilde{J}(-k)\varphi(k)\Big{]},$

where $k^{2}=|\textbf{k}|^{2}-\omega^{2}$, and $\tilde{J}(k)$ is the Fourier
transform of $J(x)$. The trick to use here for calculating the path integral
(103) with the action above is the same as the one we employed in Sect. III.1,
namely we need to discretise the integral appearing in the action $S_{0}$, and
consider a discrete set of field configurations, such that
$\int\,\mathcal{D}\phi$ can be written as a product of integrals. If we do so,
we can then use the fact that the integrand above is Gaussian in $\varphi(k)$
and repeatedly apply Gaussian integration to it to get the final result.

Alternatively, we can operate a change of path integration, by performing the
change of variables $\chi(k)=\varphi(k)-\tilde{J}(k)/(k^{2}+m^{2})$, such that
$\mathcal{D}\chi=\mathcal{D}\phi$, to isolate, in the action above, a term
that only depends on the field $\chi$, from a term independent from it, _i.e._
,

$\displaystyle S_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int\,d^{4}k\Big{[}-\chi(k)(k^{2}+m^{2}-i\varepsilon)\chi(-k)\Big{]}$
(108) $\displaystyle+$
$\displaystyle\frac{1}{2}\int\,d^{4}k\,\frac{\tilde{J}(k)\tilde{J}(-k)}{k^{2}+m^{2}}.$

If we now substitute the expression above into Eq. (103), together with
$\mathcal{D}\chi=\mathcal{D}\phi$ from the change of field variables, we
obtain the following result

$\displaystyle Z_{0}(J)$ $\displaystyle=$
$\displaystyle\exp{\left[\frac{i}{2}\,\int\,\frac{d^{4}k}{(2\pi)^{4}}\,\frac{\tilde{J}(k)\tilde{J}(-k)}{k^{2}+m^{2}}\right]}\int\,\mathcal{D}\chi\,\exp{\Big{\\{}}$
(109) $\displaystyle-$ $\displaystyle
i\int\frac{d^{4}k}{(2\pi)^{4}}\left[\chi(-k)\frac{k^{2}+m^{2}}{2}\chi(k)\right]\Big{\\}}.$

Notice, how the first term has been brought out of the path integral, as it is
independent on $\chi$. Notice moreover, that the path integral can now be
evaluated, using a generalisation of the Gaussian integration formula (see
Appendix A), and it essentially amounts to a normalisation factor (which, for
convenience, will be henceforth neglected). We then arrive at the final result
for the path integral, which can be written in a very similar form to that of
a harmonic oscillator, as

$\displaystyle Z_{0}(J)$ $\displaystyle=$
$\displaystyle\exp{\left[\frac{i}{2}\,\int\,\frac{d^{4}k}{(2\pi)^{4}}\,\frac{\tilde{J}(k)\tilde{J}(-k)}{k^{2}+m^{2}-i\epsilon}\right]}$
(110) $\displaystyle=$
$\displaystyle\exp{\left[\frac{i}{2}\,\int\,d^{4}x\,d^{4}y\,J(x)G(x-y)J(y)\right]},$

where

$G(x-y)=\int\,\frac{d^{4}k}{(2\pi)^{4}}\,\frac{e^{ik(x-y)}}{k^{2}+m^{2}-i\epsilon},$
(111)

and a small imaginary part $\varepsilon\ll 1$ has been introduced to
regularise the $\omega$-integral and allow its calculation as a contour
integral in the complex plane Byron and Fuller (1992). The quantity $G(x-y)$
is the propagator (Green’s function) for the free scalar field, which solves
the correspondent equation of motion (derived from the Euler-Lagrange
equations of the field), _i.e._ ,

$\left(\nabla^{2}-\partial_{t}^{2}+m^{2}\right)G(x-y)=\delta(x-y).$ (112)

The Green’s function can also be interpreted as the vacuum-to-vacuum two-point
correlation function between the field at point $x$ and the field at point $y$
as follows Srednicki (2007)

$G(x-y)=i\langle 0|\hat{T}\phi(y)\phi(x)|0\rangle,$ (113)

where $\hat{T}$ is the usual time-ordering operator. The two-point correlation
function is defined in terms of path integrals as follows

$\langle
0|\hat{T}\phi(x)\phi(y)|0\rangle=\int\,\mathcal{D}\phi\,\phi(x)\,\phi(y)\,\exp{\left[iS_{0}(\phi,J)\right]}.$
(114)

From this definition, by noticing that, in analogy with the derivative of the
exponential of a function we can write Kleinert (2007)

$\phi(y)e^{i\,S_{0}(\phi,J)}=-i\frac{\delta}{\delta
J(y)}e^{i\,S_{0}(\phi,J)},$ (115)

where the symbol $\delta/\delta J(y)$ stands for the functional derivative,
and it is a generalisation of the concept of derivative to functionals
Kleinert (2007) defined by the operative relation

$\frac{\delta f(x)}{\delta f(y)}=\delta(x-y).$ (116)

With this in mind, then, we can formally introduce in Eq. (113) the
substitution $\phi(x)\rightarrow-i\delta/\delta J(x)$, thus obtaining an
alternative definition for the propagator of a free field as

$\displaystyle\langle 0|\hat{T}\phi(y)\phi(x)|0\rangle$ $\displaystyle=$