Zeta function regularization, heat kernels, and one-loop determinants

In the path integral formulation of quantum field theory, a free theory is one whose Euclidean action can be expressed as
$$S[\phi] = \frac{1}{2} \int_{\mathcal{M}} \boldsymbol{\epsilon} \phi L \phi, \tag{1}$$
where $\phi$ is some tensor field and $L$ is a self-adjoint differential operator whose spectrum is bounded below. For example, if we make $\phi$ a real scalar field and $L$ the differential operator
$$L = m^2 - \nabla^2,$$
then $S[\phi]$ is the action of a free scalar with mass $m.$ It is useful to have general techniques for studying the path integrals of free field theories, not only because free theories serve as an important testing ground for interesting physical ideas, but also because in the semiclassical limit, any quantum field theory can be approximated as free. In particular, if we consider some Euclidean quantum field theory defined by the path integral
$$\int \mathcal{D} \phi e^{- \frac{1}{\hbar} S[\phi]},$$
then in the limit $\hbar \rightarrow 0$ the path integral can be approximated as a sum over stationary points of the action — these are the classical configurations of the theory — together with a perturbative-in-$\hbar$ series of quantum corrections, the path integral of fluctuations about the classical configurations. The first nontrivial term in this expansion will be quadratic in the field perturbations, since the linear term must vanish due to the assumption that the action is stationary. By applying integration by parts with appropriate boundary conditions, a quadratic functional of the fields can always be put in the form (1), with $\phi$ replaced by a linearized perturbation $\delta \phi$.

This note explains the notion of a one-loop determinant, which is a way of thinking of the path integral of a free theory as computing the (renormalized) determinant of a differential operator. It discusses how to renormalize this determinant using zeta function regularization, and how to compute that regularized determinant in the abstract using a heat kernel.

In section 1, I show how the path integral of a free theory can be interpreted as computing the determinant of a differential operator, which explains the terminology "one-loop determinant." (The term "one-loop" comes from the fact that in the semiclassical limit, the "free theory approximation" shows up at the first nontrivial order in $\hbar$.)

In section 2, I introduce zeta function regularization and show how it is related to the heat kernel $K(t)$ of $L.$

In section 3, I use zeta function regularization to compute the partition function of the quantum mechanical harmonic oscillator.

In section 4, I give some mathematical comments about the asymptotics of the heat kernel, and explain what I think are the failings of the (sort-of-incorrect) formula that is usually written down,
$$Z = \exp\left[ \frac{1}{2} \int_{0}^{\infty} dt\, \frac{1}{t} K(t) \right].$$
I also comment briefly on the "physicalness" of zeta function regularization by relating it to other regularization techniques, and refer to a blog post of Terry Tao for further information.

I learned this material from slightly unconventional sources. I learned about zeta function regularization and heat kernels from Stephen Hawking's essay in this "Einstein centenary" collection. I think I learned how to apply zeta function regularization to the harmonic oscillator from this paper, but I'm not sure; I learned this material >2 years ago and am adapting this post from notes I took at the time. While fixing some details in this post I read Hawking's original paper on the subject, which is fantastic. I found section 2 of this paper helpful for understanding how to think about Hawking's infamous number $\mu$ that appears in the path integral measure. I also highly recommend Mark Kac's delightful paper "Can one hear the shape of a drum?", which gives a great, intuitive explanation of the asymptotic properties of the heat kernel.

Prerequisites: Path integral formulation of quantum field theory. The discussion of asymptotic properties of the heat kernel in section 4 uses abstract index notation, and requires knowing enough differential geometry to understand covariant derivatives. You don't need to know differential forms to understand this post, but I'll always write $\int_{\mathcal{M}} \boldsymbol{\epsilon}$ instead of $\int d^d x \sqrt{g}.$

Table of Contents

1. Path integrals and determinants
2. Zeta function regularization
3. The harmonic oscillator partition function
4. Some mathematical comments
1. Asymptotics of the heat kernel
2. Smoothed sums
   

1. Path integrals and determinants

Suppose we have a Euclidean field theory defined by the action $S[\phi]$ on the manifold $\mathcal{M}$ with metric $g_{ab}.$ For simplicity, we'll assume the fields are real, but there's no major conceptual difference in working with theories of complex fields. The action is assumed to be free, meaning it can be written in the form
$$S[\phi] = \frac{1}{2} \int_{\mathcal{M}} \boldsymbol{\epsilon} \phi L \phi$$
for some self-adjoint differential operator $L$ whose spectrum is bounded below. We want to evaluate the partition function of the theory, which is the path integral
$$Z = \int \mathcal{D}\phi e^{- S[\phi]}.$$
The trick will be to use the fact that $L$ is self-adjoint to find an orthonormal eigenbasis of field configurations, i.e., a set of field configurations $\{\phi_n\}$ satisfying
$$L \phi_n = \lambda_n \phi_n$$
and
$$\int_{\mathcal{M}} \boldsymbol{\epsilon} \phi_n \phi_m = \delta_{n m}.$$
If we expand a generic field configuration as
$$\phi = \sum_n c_n \phi_n,$$
then it seems natural to think of the path integral measure $\mathcal{D}\phi$ as being expressible as
$$\mathcal{D}\phi = \prod_{n} \int d c_n.$$
In fact, there are times when interpreting the path integral measure this way gives inconsistent physical answers. Many of these are ameliorated by instead writing
$$\mathcal{D}\phi = \prod_{n} \frac{1}{\sqrt{2 \pi \mu}} \int d c_n,$$
with $\mu$ an as-yet unspecified constant, and fixing $\mu$ at the end of the calculation by imposing an appropriate consistency condition. I don't think it should be so surprising that such a modification is occasionally necessary. There's no real reason to think the path integral measure, this unwieldy behemoth that we don't understand on rigorous footing, should be so simple as $\mathcal{D} \phi = \prod_n \int d c_n.$ We should allow ourselves some freedom in that identification; adding a constant $\mu$ to each mode integral is a simple way to do that which is "ultralocal" in frequency space, and it turns out to be enough in many cases to get useful and interesting physical answers. 

If we expand the action in terms of the eigenbasis $\{\phi_n\},$ we find
$$S[\phi] = \frac{1}{2} \sum_{m, n} c_m c_n \lambda_n \int_{\mathcal{M}} \boldsymbol{\epsilon} \phi_m \phi_n = \frac{1}{2} \sum_{n} \lambda_n c_n^2.$$
The path integral
$$\int \mathcal{D}\phi e^{- S[\phi]}$$
then turns into a product of Gaussian integrals for all nonzero eigenvalues, plus some overall prefactor corresponding to the zero eigenvalues. We have
$$\int \mathcal{D}\phi e^{- S[\phi]} = \left(\frac{1}{\sqrt{2 \pi \mu}}\right)^{\text{# zero modes}} \operatorname{Vol}(\text{zero modes}) \prod_{n} \frac{1}{\sqrt{2 \pi \mu}} \int d c_n e^{- \frac{1}{2} \lambda_n c_n^2}. \tag{2}$$
Here $\operatorname{Vol}(\text{zero modes})$ means the volume of the domains of the coefficients of the "zero modes," i.e., the eigenfunctions $\phi_n$ with eigenvalue zero. There may be no such eigenfunctions, in which case this factor is set to $1.$ There may be zero modes whose coefficients take arbitrary values, in which case the corresponding prefactor is thought of as an unphysical divergence that should be renormalized away. In some cases, though, the prefactor $\operatorname{Vol}(\text{zero modes})$ is a finite quantity that depends on the geometry of spacetime in a nontrivial way; in these cases, the zero-mode prefactor makes an essential contribution to the path integral. We won't see an example of that in this post, but maybe in a future one; these zero modes play an essential role, for example, in establishing higher-genus modular invariance of the two-dimensional compact boson CFT.

By the way, generally speaking the prefactor $\sqrt{2 \pi \mu}^{- \text{# zero modes}}$ is physically unimportant. If the number $\mu$ doesn't depend on the geometry or on the parameters of the theory, and if the number of zero modes is a topological invariant, which is often the case, then $\sqrt{2 \pi \mu}^{-\text{# zero modes}}$ has no dependence on the parameters of the theory and therefore has no physical meaning; it's like a constant shift of the ground state energy of the theory.

The Gaussian integrals in (2) can be evaluated in the standard way, giving the result
$$\int \mathcal{D} \phi e^{- S[\phi]} = \left(\frac{1}{\sqrt{2 \pi \mu}}\right)^{\text{# zero modes}} \operatorname{Vol}(\text{zero modes}) \prod_{n} \frac{1}{\sqrt{\mu \lambda_n}}.$$
Here the product ranges over all nonzero eigenvalues of $L.$ Formally, this can be rewritten as
$$\left(\frac{1}{\sqrt{2 \pi \mu}}\right)^{\text{# zero modes}} \operatorname{Vol}(\text{zero modes}) \det \left(\mu L_{\times} \right)^{-1/2},$$
where the notation $L_{\times}$ refers to the restriction of $L$ to its nonzero eigenvectors. We will drop this notation going forward, though, and just ask the reader to keep in mind that whenever we take a determinant of $L$, or write any function of its eigenvalues, we really mean only the nonzero eigenvalues. 

The differential operator $\mu L_{\times}$ is unbounded, so it doesn't actually have a determinant. But there's a convenient and universal way to assign a finite value to this determinant, which seems to have useful physical consequences. This is the subject of the next section.

One comment, though, before we proceed: everything I've written so far, and everything I'll write below, assumes that the eigenvalue spectrum of $L$ is discrete. This is a feature of compact manifolds. For non-compact manifolds, the eigenvalue spectrum is often continuous, but also the assumption made above that a generic field configuration can be written as $\phi(x) = \sum_{n} c_n \phi_n(x)$ for arbitrary values of $c_n$ is false; boundary conditions on the field, like requiring it to fall off at a certain rate near infinity, constrain the $c_n$ coefficients relative to one another. This is one of many "infrared" issues that appear in trying to study quantum field theory on non-compact manifolds. One hopes, however, that non-compact manifolds of interest can be achieved as limits of compact manifolds, and that the properties of a quantum theory on a non-compact manifold can be computed as limits of properties of that theory on compact manifolds; this is what happens, for example, when you study quantum field theory in Minkowski spacetime by defining it on a spacelike sphere or torus and taking the limit as the size of that sphere/torus goes to infinity.

2. Zeta function regularization

We will now start assuming that $L$ is a positive semidefinite operator, i.e., that all of its nonzero eigenvalues are positive. This is not an especially significant restriction, as any operator whose spectrum is bounded below can be made positive semidefinite by a constant shift.

Our goal is to assign a finite value to the formal expression
$$\det\left( \mu L \right) = \prod_{n} \mu \lambda_n.$$
If this product were convergent, then we could take its logarithm and obtain
$$\log\det\left(\mu L \right) = \sum_n \log\left(\mu \lambda_n \right). \tag{3}$$
A convenient way to rewrite the logarithm of a number is as the linear contribution to its power function, i.e., as
$$\log(x) = - \left. \frac{d}{ds} x^{-s} \right|_{s=0}.$$
Substituting this into equation (3) gives
$$\log\det\left(\mu L \right) = - \sum_n \frac{d}{ds} \left. (\mu \lambda_n)^{-s} \right|_{s=0}.$$
If this series were uniformly convergent, then we could exchange the sum and the derivative to obtain
$$\log\det\left(\mu L \right) = - \left. \frac{d}{ds} \left( \sum_n (\mu \lambda_n)^{-s} \right)\right|_{s=0}.$$
Consider now $s$ as a complex number. For values of $s$ with sufficiently large real part, unless the eigenvalues $\{\lambda_n\}$ are particularly badly behaved at large values of $n,$ the sum
$$\sum_n (\mu \lambda_n)^{-s}$$
converges. In fact, it converges uniformly on compact subsets in the $s$-plane, and therefore by Morera's theorem is an analytic function on its domain of definition. On this domain, we define the $\zeta$-function
$$\zeta_{\mu L}(s) = \sum_{n} (\mu \lambda_n)^{-s}.$$
It can be extended to other parts of the complex plane by analytic continuation. If this extension is unique and regular at $s=0$, then we can formally define the determinant of $\mu L$ by
$$\det\left(\mu L \right) = e^{- \zeta_{\mu L}'(0)}.$$

Supposedly, one can show that when $L$ is a so-called elliptic differential operator, the asymptotic growth of the eigenvalues for large $n$ is such that the $\zeta$ series is guaranteed to converge uniformly for $\operatorname{Re}(s)$ above some threshold, and such that the analytic continuation is guaranteed to be regular at $s=0.$ For this result, Hawking cites the essay "Complex powers of an elliptic operator" in this book, but the paper is pretty technical and I haven't been able to understand it. No big deal, though; we can check that the eigenvalues of any particular operator we're interested in admit an appropriate analytic continuation, and indeed we'll do this for the harmonic oscillator in section 3. In fact, we'll see in section 4 that the question of whether the $\zeta$ series converges for large enough $\operatorname{Re}(s)$ can be framed in terms of the asymptotic behavior of a particular function of the eigenvalues, and give an intuitive explanation for why we should expect it to converge.

In the zeta function's domain of definition, we can write
$$\zeta_{\mu L}(s) = \mu^{-s} \sum_{n} \lambda_n^{-s} = \mu^{-s} \zeta_L(s).$$
Because $\mu^{-s}$ is regular everywhere in the $s$ plane, this expression holds everywhere. So in particular, near $s=0,$ we have
$$\zeta_{\mu L}(s) = (1 - s \log\mu) \zeta_{L}(s),$$
and hence
$$\zeta_{\mu L}'(0) = - \zeta_{L}(0) \log\mu + \zeta_{L}'(0).$$
It is often the case, though not always, that $\zeta_{L}(0)$ vanishes; when this happens, we have $\zeta_{\mu L}'(0) = \zeta_{L}'(0),$ and the normalization constant $\mu$ makes no physical contribution to the path integral. In any case, it is clear that we only need to study the function $\zeta_L(s)$ in a neighborhood of $s=0,$ since $\zeta_{\mu L}'(0)$ is completely determined by $\zeta_{L}(0),$ $\zeta_{L}'(0),$ and $\mu.$

There is a useful way to rewrite the function $\zeta_{L}(s)$, which is to insert the gamma function
$$\Gamma(s) = \int_{0}^{\infty} dt\, t^{s-1} e^{-t}.$$
We can write
$$\zeta_{L}(s) = \frac{1}{\Gamma(s)} \sum_{n} \int_{0}^{\infty} dt\, t^{s-1} e^{-t} \lambda_n^{-s}.$$
Making the substitution $t \mapsto t \lambda_n$ gives
$$\zeta_{L}(s) = \frac{1}{\Gamma(s)} \sum_{n} \int_{0}^{\infty} dt\, t^{s-1} e^{-t \lambda_n}.$$
For well-behaved $\lambda_n$ and sufficiently large values of $\operatorname{Re}(s),$ the sum $\sum_n t^{s-1} e^{-t \lambda_n}$ converges uniformly on compact subsets of $t$, which means we can interchange the sum and the integral to write
$$\zeta_L(s) = \frac{1}{\Gamma(s)} \int_{0}^{\infty} dt\, t^{s-1} \sum_n e^{-t \lambda_n}.$$
(The only real restriction to guarantee this uniform convergence is that $\operatorname{Re}(s)$ needs to be large enough to cancel any small-$t$ poles in $\sum_n e^{-t \lambda_n}$; more on that in section 4.)

The sum $\sum_n e^{-t \lambda_n}$ can in principle be evaluated directly if we know all of the eigenvalues of the operator $L,$ which will then let us evaluate $\zeta_L(s)$ near $s=0.$ Indeed, we'll do this explicitly in section 3 for the harmonic oscillator partition function. But it isn't always the case that we know all of the eigenvalues, or how to sum them, or how to understand the integral of $t^{s-1}$ against that sum. But we don't actually need all of the eigenvalues to evaluate $\zeta_L(s)$; we only need the sum. The total sum has an elegant relationship to the heat kernel of the operator $L,$ which in some cases can be studied using techniques from partial differential equations even if we can't study the eigenvalues directly.

The idea is to write the sum $\sum_n e^{- t \lambda_n}$ in terms of an orthonormal basis of eigenfunctions $\{\phi_n\}$ as
$$\sum_n e^{-t \lambda_n} = \int_{\mathcal{M}, x} \boldsymbol{\epsilon} \sum_n e^{- t \lambda_n} \phi_n(x) \phi_n(x).$$
If we define the function $K(x, x', t)$ by
$$K(x, x', t) = \sum_n e^{- t \lambda_n} \phi_n(x) \phi_n(x'),$$
then we can express this as
$$\sum_n e^{- t \lambda_n} \int_{\mathcal{M}, x} \boldsymbol{\epsilon} K(x, x, t).$$
The function $K(x, x', t)$ for $t > 0$ is interesting because it solves the heat equation
$$\frac{\partial K(x, x', t)}{\partial t} + L_{x} K(x, x', t) = 0,$$
which is straightforward to check. The heat equation is linear and involves only a single derivative in $t,$ so a solution of the heat equation is determined uniquely for all $t > 0$ by its boundary condition at $t=0.$ In the limit $t \rightarrow 0,$ the function $K(x, x', t)$ limits to
$$K(x, x', 0) = \sum_n \phi_n(x) \phi_n(x') = \delta(x, x'),$$
where the delta function is normalized with respect to the metric, so that it satisfies
$$\int_{\mathcal{M}, x} \boldsymbol{\epsilon} \delta(x, x') f(x) = f(x').$$
Because $K(x, x', t)$ is the unique solution to the heat equation for the differential operator $L$ with boundary condition $K(x, x', 0) = \delta(x, x'),$ it is called the heat kernel for $L.

Because of the connection to the heat kernel, we'll now start using the notation
$$K(t) = \sum_n e^{- t \lambda_n}.$$
This is the "trace of the heat kernel," $K(t) = \int_{\mathcal{M}, x} \boldsymbol{\epsilon} K(x, x, t).$ Per our previous discussion, the zeta function for the operator $L$ is given by
$$\zeta_{L}(s) = \frac{1}{\Gamma(s)} \int_{0}^{\infty} dt\, t^{s-1} K(t). \tag{4}$$

3. The harmonic oscillator partition function

The harmonic oscillator can be thought of as a one-dimensional quantum field theory, where the one-dimensional base space $\mathcal{M}$ encodes the time variable of the theory and the field is the position function $x : \mathcal{M} \rightarrow \mathbb{R}.$ In Lorentz signature, the harmonic oscillator is studied on the base space $\mathcal{M} = \mathbb{R},$ with action
$$S = \int dt\, \left( \frac{m}{2} x'(t)^2 - \frac{m \omega^2}{2} x(t)^2 \right).$$
In Euclidean signature, the action is
$$S = \int dt\, \left( \frac{m}{2} x'(t)^2 + \frac{m \omega^2}{2} x(t)^2 \right).$$
The Euclidean theory is the right one to use to compute the partition function of the harmonic oscillator; because what we want to compute is
$$Z(\beta) = \operatorname{tr}(e^{- \beta H}),$$
we can compute this as the Euclidean time evolution operator on a periodic manifold of length $\beta.$ (Length $\beta$ because we want to evolve for Euclidean time $\beta$, periodic because we want the path integral to implement the trace.) In other words, we want to perform the path integral of the harmonic oscillator over the circle $S^1$ with metric
$$ds^2 = \left(\frac{\beta}{2 \pi} \right)^2 d\theta^2.$$

After integrating by parts, the Euclidean action for the harmonic oscillator on this manifold can be written as
$$S[x] =\frac{1}{2} \int_{0}^{\beta} dt\, x(t) \left[ m \omega^2 - m \partial_t^2 \right] x(t).$$
This is a free theory with differential operator $L = m \omega^2 - m \partial_t^2.$ On a circle of length $\beta,$ the normalized eigenfunctions of this differential operator are the Fourier modes
$$x_{n}(t) = \frac{1}{\sqrt{\beta}} e^{2 \pi i n t/\beta}$$
with eigenvalues
$$\lambda_n = m \left( \omega^2 + \frac{4 \pi^2 n^2}{\beta^2} \right).$$
The trace of the heat kernel for this differential operator is thus
$$K(t) = e^{-m t \omega^2} \sum_{n=-\infty}^{\infty} \exp\left[ - \frac{4 \pi^2 m t}{\beta^2}n^2 \right].$$
The sum in this expression is a very famous function called the Jacobi theta function. We can study it by applying the Poisson summation formula, which gives us the expression
$$K(t) = e^{- m t \omega^2} \frac{\beta}{2 \sqrt{\pi m t}} \sum_{n=-\infty}^{\infty} \exp\left[ - \frac{\beta^2}{4 m t} n^2 \right].$$
The advantage of this expression is that the integral
$$\Gamma(s) \zeta_{L}(s) = \int_0^{\infty} dt\, t^{s-1} K(t) = \frac{\beta}{2 \sqrt{\pi m}} \sum_{n=-\infty}^{\infty} \int_0^{\infty} dt\, t^{s-3/2} \exp\left[ - m t \omega^2 - \frac{\beta^2}{4 m t} n^2 \right]$$
is regular even in the limit $s \rightarrow 0$ (except for the $n=0$ term, which we'll treat separately). We could ask Mathematica to do the integral for us, but it's kind of a fun exercise, so let's do it ourselves. It'll be convenient to be able to take $n$ to be positive in the following, so we'll rewrite our sum as
$$\Gamma(s) \zeta_{L}(s) = \frac{\beta}{2 \sqrt{\pi m}} \left( 2 \sum_{n=1}^{\infty} \int_0^{\infty} dt\, t^{s-3/2} \exp\left[ - m t \omega^2 - \frac{\beta^2}{4 m t} n^2 \right] + \int_{0}^{\infty} dt\, t^{s-3/2} e^{- m t \omega^2} \right).$$
The second integral can be evaluated directly as
$$\int_{0}^{\infty} dt\, t^{s-3/2} e^{- m t \omega^2} = \Gamma\left(s - \frac{1}{2} \right) (m \omega^2)^{\frac{1}{2} - s}.$$
This is analytic at $s=0,$ and evaluates to $- 2 \omega \sqrt{\pi m}.$

The remaining integral we want to evaluate is
$$I = \int_{0}^{\infty} dt\, \frac{1}{t^{3/2}} \exp\left[ - m \omega^2 t - \frac{\beta^2 n^2}{4 m} \frac{1}{t} \right].$$
To start, let's try to get this integral to look a little more Gaussian by substituting $t \mapsto u^2,$ obtaining
$$I = 2 \int_{0}^{\infty} du\, \frac{1}{u^2} \exp\left[ - m \omega^2 u^2 - \frac{\beta^2 n^2}{4 m} \frac{1}{u^2} \right].$$
The exponent is a bit alluring; because it includes factors of both $u^2$ and $u^{-2},$ there is some number $\alpha$ for which the substitution $u \mapsto \alpha/u$ will keep the exponent fixed. This value can easily be checked to be $\alpha = \beta n / 2 m \omega.$ Under this substitution, we obtain
$$I = \frac{4 m \omega}{\beta n} \int_{0}^{\infty} du\, \exp\left[ - m \omega^2 u^2 - \frac{\beta^2 n^2}{4 m} \frac{1}{u^2} \right].$$
(Here for the first time we have used the assumption $n > 0,$ which keeps the bounds of integration in the right order.) Summing the two expressions for $I$, we obtain
$$I = \frac{2 m \omega}{\beta n} \int_{0}^{\infty} du\, \left( 1 + \frac{\beta n}{2 m \omega} \frac{1}{u^2} \right) \exp\left[ - m \omega^2 u^2 - \frac{\beta^2 n^2}{4 m} \frac{1}{u^2} \right].$$
This encourages us to try substituting in a variable $w$ satisfying
$$dw = du\, \left(1 + \frac{\beta n}{2 m \omega} \frac{1}{u^2} \right).$$
This is accomplished by the change of variables
$$w = u - \frac{\beta n}{2 m \omega} \frac{1}{u}.$$
Because we have assumed $n>0,$ this variable ranges from $-\infty$ at $u=0$ to $+\infty$ at $u=+\infty.$ Furthermore, one can easily check the relation
$$m \omega^2 w^2 = m \omega^2 u^2 + \frac{\beta^2 n^2}{4 m} \frac{1}{u^2}- \beta n \omega,$$
which is exactly the exponent in our most recent integral expression for $I$ up to a constant term. So, substituting $w$ for $u$ we obtain
$$I = \frac{2 m \omega}{\beta n} \int_{-\infty}^{\infty} dw\, \exp[- m \omega^2 w^2 - \beta n \omega].$$
This is a Gaussian integral, and it evaluates to
$$I = e^{- \beta n \omega} \frac{2 \sqrt{\pi m}}{\beta n}.$$
This gives us the $s=0$ limit of the zeta function as
$$\lim_{s \rightarrow 0} \Gamma(s) \zeta_{L}(s) = \frac{\beta}{2 \sqrt{\pi m}} \left( 2 \sum_{n=1}^{\infty} e^{- \beta n \omega} \frac{2 \sqrt{\pi m}}{\beta n} - 2 \omega \sqrt{\pi m} \right) = 2 \sum_{n=1}^{\infty} \frac{e^{- \beta n \omega}}{n} - \beta \omega.$$

The sum can be evaluated by taking the derivative with respect to $\beta \omega,$ which results in an ordinary geometric series that can be summed with the usual formula. The result can then be integrated with respect to $\beta \omega,$ with the constant of integration fixed by the condition that the integral should vanish in the limit $\beta \omega \rightarrow \infty.$ This gives the expression
$$\lim_{s \rightarrow 0} \Gamma(s) \zeta_L(s) = -\frac{1}{\beta} \log(e^{\beta \omega} - 1).$$
In the small-$s$ limit, the gamma function has a pole $\Gamma(s) \sim 1/s.$ So because the right-hand side of this expression is regular, we must have $\zeta_L(0) = 0.$ As such, the renormalization constant $\mu$ discussed in the previous section makes no contribution to the one-loop determinant, and we have
$$\zeta_{\mu L}'(0) = \zeta_{L}'(0).$$
Again using $\Gamma(s) \sim 1/s,$ we obtain the identity
$$\zeta_L'(0) = \beta \omega - 2 \log(e^{\beta \omega} - 1).$$
From the formulas of the previous section, we see that the partition function of the harmonic oscillator is
$$Z[\beta] = \left(\frac{1}{\sqrt{2 \pi \mu}}\right)^{\text{# zero modes}} \operatorname{Vol}(\text{zero modes}) e^{\zeta_{L}'(0)/2}.$$
For nonzero frequencies $\omega \neq 0,$ the differential operator $m \omega^2 - m \partial_{t}^2$ has no zero modes, so we are left with
$$Z[\beta] = e^{\zeta_{L}'(0)/2} = \frac{1}{2 \sinh\left(\frac{\beta \omega}{2} \right)}. \tag{5}$$

Let's check our work. From the Hamiltonian formalism, we know that the energy levels of the harmonic oscillator are $\omega(n + 1/2)$ for $n=0$ to $n=\infty,$ so the partition function should be
$$Z[\beta] = e^{- \beta \omega/2} \sum_{n=0}^{\infty} e^{- \beta \omega n}.$$
This is a geometric series. The result, fully simplified, matches equation (5).

4. Some mathematical comments

This section is a bit of a grab-bag. In the first subsection, I'll comment on the asymptotics of the heat kernel, and on the "sort-of-correct" formula
$$Z = \exp\left[\frac{1}{2} \int_{0}^{\infty} dt\, \frac{1}{t} K(t) \right]$$
that is often written in quantum field theory texts. In the second subsection, I'll comment on the fact that under certain assumptions the $\zeta$-function regularization technique produces the same result as more "physical-seeming" regularization schemes.

4.1 Asymptotics of the heat kernel

Recall from sections 1 and 2 that the path integral of a free theory can be expressed in terms of $\zeta_{L}(0)$ and $\zeta_{L}'(0),$ and that $\zeta_{L}(s)$ can be written in terms of the heat kernel as
$$\zeta_{L}(s) = \frac{1}{\Gamma(s)} \int_0^{\infty} dt\, t^{s-1} K(t),$$
so long as $\operatorname{Re}(s)$ is large enough that the integral converges. Because $1/\Gamma(s)$ behaves like $s$ in the limit $s \rightarrow 0,$ we are tempted to just take the $s\rightarrow 0$ limit of the integral and write
$$\zeta_{L}'(0) = \int_{0}^{\infty} dt\, \frac{1}{t} K(t). \tag{6}$$

There are a few problems with this expression. First, applying the same logic without taking an $s$-derivative would lead us to believe that $\zeta_{L}(0)$ vanishes; while this is often the case, as in the harmonic oscillator example of the previous section, it is not always true. Second, the integral $\int_{0}^{\infty} dt\, \frac{1}{t} K(t)$ will usually not converge! If the eigenvalues are well behaved at infinity then we might expect the upper limit of the integral to converge, because $K(t)$ should decay exponentially in $t,$ but near $t\sim 0$ the integral will only converge if $K(0)$ vanishes. This is almost never the case; in fact, we'll argue below on quite general grounds that it will have a pole at small $t$ whose order is controlled by the dimension of the manifold $\mathcal{M}$ and by the order of the differential operator $L.$

So if we take equation (6) as our starting point, we need to evaluate it by introducing some cutoff and writing
$$\zeta_{L}'(0) = \int_{\epsilon}^{\infty} dt\, \frac{1}{t} K(t),$$
then introducing some counterterms to eliminate the divergences in $\epsilon.$ I don't really like this approach; that cutoff is hard for me to make sense of physically, and in reducing ourselves to this expression we completely miss the fact that some divergences are physical and some aren't. For example, the partition function of the harmonic oscillator has no physical divergences; the zeta function regularization we undertook in the previous section gave us finite expressions the whole way through. But the heat kernel for that problem still diverges as $1/\sqrt{t}$ in the small-$t$ limit, so equation (6) will have divergences that must be renormalized away even for the harmonic oscillator. On the other hand, there are divergences that have a more physical meaning; for example, the path integrals that compute entanglement entropy of spatial subregions of free quantum field theories do have physical divergences that represent large amounts of entanglement in the ultraviolet modes of the theory. These divergences are visible in the appropriate zeta function, and can be studied qualitatively by introducing a short-length cutoff in a neighborhood of the entangling surface and asking how the zeta function near $s=0$ depends on the size of the cutoff. If you start your analysis with equation (6), then I think it becomes harder to keep track of which divergences are mathematical trivialities and which ones really tell you something physical about the theory.

But that's enough ideology. You can, of course, start with equation (6) and get perfectly good physical answers by being careful about renormalization. Many people do, successfully! And I'm sure some people would argue that the UV divergence in the entanglement entropy isn't a physical property of the theory, since it is expected to smooth out when the theory is coupled to gravity. Let's drop the matter and move onward to study the small-$t$ asymptotics of the heat kernel.

We will restrict our attention to a class of differential operators that might initially seem fairly narrowly defined, but that encompasses a large number of physically interesting theories. A general $n$-th order differential operator can be expressed in terms of the covariant derivative $\nabla_a$ as
$$L \phi = \alpha_{(0)} \phi + \alpha_{(1)}^{a} \nabla_a \phi + \alpha_{(2)}^{a b} \nabla_{a} \nabla_{b} \phi + \dots + \alpha_{(n)}^{a_1 \dots a_n} \nabla_{a_1} \dots \nabla_{a_n} \phi.$$
We will restrict our attention to cases where all of the coefficients $\alpha_{(j)}^{a_1 \dots a_j}$ are tensors locally constructed from the metric; in particular, the only tensors we will consider are combinations of the metric $g_{ab},$ the Riemann tensor $R_{abc}{}^{d},$ and the $D$-dimensional volume form $\epsilon_{a_1 \dots a_D}.$ These various combinations may be multiplied by real coefficients, but not by arbitrary functions of the base manifold.

The basic idea behind understanding the heat kernel of such an operator is that physically, the heat equation for $L$ is supposed to govern diffusion on a manifold $\mathcal{M}$ where the local rate of diffusion is set by the differential operator $L.$ The heat kernel $K(x, x', t)$ is supposed to represent the $x$-profile of diffusing matter that, at time $t=0,$ was concentrated entirely at the point $x'.$ Since diffusion is local, at very small values of $t$ the heat kernel should be basically unaware of the curvature of the manifold on which $L$ is defined. So if we choose locally orthonormal coordinates at the point $x$ — so that the Christoffel symbols vanish at that point — then for small $t$ the heat kernel shouldn't know that the Christoffel symbols are nonzero away from $x,$ and should approach that of the differential operator
$$L^{(\text{flat})} \phi = \alpha_{(0)} \phi + \alpha_{(1)}^{\mu} \partial_{\mu} \phi + \alpha_{(2)}^{\mu \nu} \partial_\mu \partial_\nu \phi + \dots + \alpha_{(n)}^{\mu_1 \dots \mu_n} \partial_{\mu_1} \dots \partial_{\mu_n} \phi, \tag{5}$$
where terms in the coefficients that depend on the Riemann tensor are set to zero. In fact, the terms depending on the volume form $\epsilon_{a_1 \dots a_D}$ must also vanish, because the volume form is antisymmetric and all derivatives commute when the curvature vanishes. Since the only remaining tensor is the flat metric $\eta_{\mu \nu},$ which has two indices, and since $L^{(\text{flat})} \phi$ must have the same number of indices as $\phi,$ equation (5) can contain only even derivatives of $\phi.$

The eigenfunctions of $L^{(\text{flat})}$ are all proportional to the plane waves $e^{i p \cdot x},$ with corresponding eigenvalues
$$\alpha_{(0)} - \alpha_{(2)}^{\mu \nu} p_{\mu} p_{\nu} + \dots + (-1)^{k} \alpha_{(2k)}^{\mu_1 \dots \mu_{2k}} p_{\mu_1} \dots p_{\mu_{2k}}.$$
Because each $\alpha$ coefficient is constructed from the metric $\eta_{\mu \nu},$ this eigenvalue must be equal to
$$\alpha_{(0)} + \alpha_{(2)} p^2 + \alpha_{(2k)} p^{2 k}$$
for some appropriate real coefficients $\alpha_{(j)}.$

The number of independent eigenfunctions at each momentum value $p^{\mu}$ is equal to the number of linearly independent degrees of freedom in the tensor $\phi$; if $\phi$ has $N$ indices, then there are $D^N$ independent eigenfunctions. For example, if $\phi$ is really a rank-(2,0) tensor $\phi^{\mu \nu},$ then the eigenfunctions with eigenvalue $p$ are
$$(\phi_{p; j_1, j_2})^{\mu \nu} (x) = e^{i p \cdot x} \delta_{j_1}{}^{\mu} \delta_{j_2}{}^{\nu}.$$
So the heat kernel can be written as
$$K^{\text{flat}} (x, x', t) = \int d^{D} p\, D^N e^{-t(\alpha_{(0)} + \alpha_{(2)} p^2 + \dots + \alpha_{(2k)} p^{2k})} \frac{e^{i p \cdot (x-x')}}{(2 \pi)^D},$$
where the factor of $1/(2 \pi)^D$ comes from the requirement that the eigenfunctions be $\delta$-function normalized.

This integral converges whenever the leading coefficient in the exponential at large values of $p,$ $- \alpha_{(2k)},$ is negative. In that case, we can even estimate the behavior of $K^{\text{flat}}(x, x, t)$ for small $t.$ We first set $x' = x$ to obtain
$$K^{\text{flat}} (x, x, t) = \frac{D^N}{(2 \pi)^{D}} \int d^{D} p\, e^{-t(\alpha_{(0)} + \alpha_{(2)} p^2 + \dots + \alpha_{(2k)} p^{2k})}.$$
We then substitute $p \mapsto p\, t^{-\frac{1}{2k}}$ to obtain
$$K^{\text{flat}} (x, x, t) = \frac{D^N}{(2 \pi)^{D} t^{\frac{D}{2k}}} \int d^{D} p\, e^{- \alpha_{(0)} t - \alpha_{(2)} p^2 t^{\frac{k-1}{k}} - \dots - \alpha_{(2k)} p^{2k}}.$$
For small $t$ the integral is obviously dominated by the last term in the exponential, which is the only one with no $t$ dependence. The integral of this term can be computed by standard methods (converting to higher-dimensional polar coordinates), but we won't actually need the expression, all we need to know is that it is finite in the limit $t\rightarrow 0.$ This tells us that for small $t$, the function $K^{\text{flat}}(x, x, t)$ is dominated by a term of order $1/t^{D/2k}.$ Since we expect on physical grounds that our actual heat kernel $K$ should approach $K^{\text{flat}}$ at small $t,$ this tells us that on our base manifold $\mathcal{M},$ the trace of the heat kernel of the differential operator $L$ goes like
$$K(t) = \int_{\mathcal{M}} \boldsymbol{\epsilon} K(x, x, t) \sim \operatorname{Vol}(\mathcal{M}) \frac{\#}{t^{\frac{D}{2k}}} + \text{higher order in $t$}.$$
If $\mathcal{M}$ is not compact, then the volume factor should be treated as an infrared divergence.

We can now return to our formula (4) for the zeta function, reproduced here:
$$\zeta_{L}(s) = \frac{1}{\Gamma(s)} \int_0^{\infty} dt\, t^{s-1} K(t) \tag{6}.$$
Because we know that for small $t$ the trace of the heat kernel goes as $t^{-D/2k},$ we know that the lower limit of the integral converges for $\operatorname{Re}(s) > \frac{D}{2k}.$ Unfortunately we haven't provided a similar argument for the convergence of the integral in the $t \rightarrow \infty$ limit. For that limit of the integral to converge, it must be the case that the large eigenvalues aren't too densely packed — we want the sum $\sum_n e^{- \lambda_n t}$ to decay exponentially at large $t,$ which will happen so long as the number of eigenvalues near the value $\lambda$ grows only polynomially as $\lambda$ becomes large. But this can actually be shown from what we already have if we apply the Hardy-Littlewood tauberian theorem. What we do is rewrite the heat kernel as
$$K(t) = \sum_{n} e^{- \lambda_n t} = \int_{0}^{\infty} dx\, e^{- x t} \sum_{n} \delta(x - \lambda_n).$$
If we define $F(x)$ to be the piecewise-constant function
$$F(x) = \{n \text{ if } \lambda_n \leq x < \lambda_{n+1} \},$$
then we have
$$d F(x) = \sum_{n} \delta(x - \lambda_n) dx,$$
and hence
$$K(t) = \int_{0}^{\infty} dx\, e^{-x t} dF(x).$$
The integral formulation of the Hardy-littlewood theorem then tells us that because $K(t)$ asymptotes for small $t$ to
$$K(t) \sim_{t \rightarrow 0} \frac{c}{t^{\frac{D}{2k}}},$$
the function $F(x)$ asymptotes for large $x$ to
$$F(x) \sim_{x \rightarrow \infty} \frac{c}{\Gamma\left(\frac{D}{2k} + 1\right)} x^{\frac{D}{2k}}.$$
Since $F(x)$ counts the total number of eigenvalues with value at or below $x,$ we learn that the asymptotic rate of growth is polynomial in $x,$ and therefore sub-exponential, as desired.

4.2 Smoothed sums

In section 2, we introduced zeta function regularization as a way of assigning a finite value to an infinite product like
$$\prod_{n=1}^{\infty} \lambda_n, \tag{7}$$
where $\lambda_n$ is a nondecreasing sequence because the spectrum of our differential operator is bounded below.  This product clearly does not converge unless the eigenvalues are bounded above by $1.$

 

But let's consider for a moment what equation (7) means. It is defined as the limit
$$\lim_{N \rightarrow \infty} \prod_{n=1}^{N} \lambda_n. \tag{8}$$

But this is only one way of defining what we mean by an infinite product. There is another way, which I will shortly argue is more physical, which is to introduce a smooth regulating function $\eta(x),$ which has the properties that (i) it is compactly supported, and (ii) it equals one in a neighborhood of $x=0.$ We then define our infinite product as
$$\lim_{N \rightarrow \infty} \prod_{n=1}^{\infty} \lambda_n^{\eta\left( \frac{n}{N} \right)}. \tag{9}$$
As $N$ gets larger, more and more of the integers $n$ have the property that $n/N$ is in the neighborhood of $x=0$ with $\eta = 1.$ So as $N$ gets larger, we fully expect the limit to be the same as the limit in equation (8). In fact, for any product that is absolutely convergent, an appropriate application of the dominated convergence theorem can be used to show that (9) and (8) give identical results. But equation (9) has a serious advantage for non-convergent products, which is that it is smooth in $N,$ and therefore can be given an asymptotic expansion in the large-$N$ limit.

So here's the trick: when we study quantum field theory, we really have in mind that we're studying some low-energy theory whose ultraviolet completion we don't understand. Renormalization is the process of "fuzzing out" whatever happens at high energies, with the assumption that it is an artifact of our ultraviolet ignorance. One approach to renormalizing a path integral is to study the asymptotic expansion of equation (9) in the large-$N$ limit, and to keep only the finite-in-$N$ term of that expansion. The terms that diverge with $N$ are contributions from UV modes of the theory; if we're lucky, they can be eliminated by adding local counterterms.

I won't be any more precise than that, but I will tell you that as I understand it, the $\zeta$-function regularized result for $\prod_{n=1}^{\infty} \lambda_n$ reproduces the finite-in-$N$ term in the large-$N$ expansion of (9) under quite general assumptions. I don't know much about the general theory, but there are a few very clear examples worked out in this blog post by Terry Tao. I'll send you there for further explorations. If you find a good, clean way to explain why this should be true in general, send me a message!