Unitarity Violation in Field Theories of Lee-Wick’s Complex Ghost

Theories with fourth-order derivatives, including the Lee-Wick ﬁnite QED model and Quadratic Gravity, have a better UV behaviour, but the presence of negative metric ghost modes endanger unitarity. Noticing that the ghost acquires a complex mass by radiative corrections, Lee and Wick, in particular, claimed that such complex ghosts would never be created by collisions of physical particles because of energy conservation, so that the physical S-matrix unitarity must hold. We investigate the unitarity problem faithfully working in the operator formalism of quantum ﬁeld theory. When complex ghosts participate, a complex delta function (generalization of Dirac delta function) appears at each interaction vertex, which enforces a speciﬁc conservation law of complex energy. Its particular property implies that the naive Feynman rule is wrong if the four-momenta are assigned to the internal lines after taking account of the conservation law in advance. We show that the complex ghosts are actually created and unitarity is violated in such fourth-order derivative theories. We also ﬁnd a deﬁnite energy threshold below which the ghosts cannot be created: The theories are unitary and renormalizable below the threshold.


I. INTRODUCTION
It has long been a challenge to construct a consistent quantum filed theory (QFT) of gravity, though a well-known candidate is Quadratic Gravity (QG) whose renormalizability was proven by Stelle [1] long ago.The action of QG contains in addition to the Einstein-Hilbert term R, the Weyl tensor squared C 2 µνρσ as well as R 2 (see e.g.[2,3] for a review).The canonical dimension of these two terms is four, which means that any renormalizable QFT of gravity will contain them [4].The main reason that QG is renormalizable is that it has a built-in ultraviolet (UV) regulator, the spin-two ghost field [1], which however endangers the unitarity of the theory.The origin of the good side and the bad side is the same; fourthorder derivatives in the equations of motion.The propagator of a fourth-order theory can be written as a sum of two propagators of a second-order theory, where one of them must come with an opposite sign, where the wrong-sign propagator originates from the wrong sign of the kinetic term.This implies that there exists a negative-metric ghost. 1 Our interest in this paper is focused on this troublesome feature, and we shall investigate the unitarity problem in great detail.
Long before Stelle's proof for renormalizability of QG, Lee and Wick [8][9][10] used the aforementioned built-in mechanism to soften the UV behavior of QFT.They considered a QED-like theory [10] with a spin-one massive ghost field and made a very remarkable observation: The radiative corrections to the ghost propagator due to light fermion (e.g., lepton pair) loops makes the single pole split into two complex conjugate poles.That is, the mass and energy of the ghost become complex.This fact became an essential part of their proof of the unitarity: Ghosts cannot be produced in a scattering process of physical particles possessing real energies because of the conservation of energy (in particular, because of the conservation of the imaginary part of the energy). 2 If the ghosts are not produced from any physical initial states, the unitarity of the S-matrix restricted to the physical particles alone, i.e., the physical unitarity, holds.Although this physical intuition sounds correct, it seems that Feynman diagrams do not share this intuition since one has to give a precise prescription how to integrate internal momenta [8][9][10]14].In fact, how to choose an appropriate integration contour in the complex energy plane has been the target of theoretical physicists since then and even recently [15][16][17][18].3 The approach described above is based on the standpoint that a theory is defined by its Feynman expansion [11].
Our interest in this paper is to investigate the aforementioned unitarity problem solely based on QFT.We will use the tools of QFT, no more, no less.How to integrate the loop momenta in a Feynman diagram, apart from its regularization, is dictated by QFT with no room left for speculation.If energy is complex, we encounter a mathematical expression, a complex distribution, which to our knowledge has not been sufficiently explored so far aside from a brief discussion in Nakanishi's work [26][27][28].It is just the integral lim T →∞ 1 2π which appears at each vertex, where E is the sum i E i of the energies of particles i entering the vertex.If the energy E is real, the limit T → ∞ is the Dirac delta function, as we all know.But what is the limit if E is complex?Nakanishi [26][27][28] called it a complex delta function and derived some of its mathematical properties, but he did not recognize that the complex distribution plays an essential role in investigating unitarity.The integral (1) with a complex E appears at each vertex of the Dyson expansion in perturbation theory if ghost lines are attached.Usually, we do not care about this integral because it just gives a delta function which expresses the energy conservation at each vertex.Thus, the mathematical property of the integral (1) is intimately related to the energy conservation, especially in the case that the energy is complex.Or, speaking more strongly, the property of the complex delta function (or distribution) defined by the integral ( 1) is all of the content of the "energy conservation law" in the case where complex ghosts participate.
We will see that δ c (E) is non-vanishing even when Im E = 0, thus implying that ghosts can be produced with finite (i.e., non-zero) probability through a collision of physical particles.This leads us to the conclusion that the physical unitarity of the S-matrix in fourth-order derivative theories such as the Lee-Wick finite QED and QG is violated in QFT. 4 A positive message, however, is that these theories are good effective theories in the sense that the physical unitarity of the S-matrix is satisfied below a definite threshold energy for the ghost production, where internal ghost propagators in a Feynman diagram may still be present. 5 In connection with the property of the complex delta function, we should mention a very important point which invalidates the usual Feynman rule with energy-momentum conservation used at each vertex in advance.For instance, consider a self-energy type oneloop diagram consisting of a physical particle ψ with real mass squared µ 2 and a ghost ϕ with complex mass squared M 2 .Many would naively write down the expression Σ(p) ∝ d 4 qD ϕ (q) D ψ (p − q) = d 4 q 1 q 2 + M 2 but this is wrong!T. D. Lee also started with this expression and correctly calculated this integral, even taking the deformed q 0 -integration contour into account properly, and reached the conclusion that the amplitude has no imaginary part so that the production of a ghost ϕ and a physical particle ψ would not occur.This is an incorrect conclusion obtained from a correct calculation but from the wrong starting expression.The correct expression is not Eq.( 2), but The k 0 -integration of the complex delta function δ c (k 0 +q 0 −p 0 ) does not give the substitution rule k 0 → p 0 − q 0 alone when the multiplied function is not analytic in k 0 .In the usual 4 By the physical unitarity we mean that the probability interpretation of quantum theory is possible.The physical unitarity of the S-matrix should be distinguished from the total S-matrix unitarity (S † S = SS † = 1) which is satisfied if the Hamiltonian is hermitian.In the following discussions we sometimes suppress "physical" if it does not lead to a confusion. 5Because of the imaginary part of the complex mass, the meaning of the threshold should be slightly modified.
Feynman graph case, the multiplied function is the propagator, here D ψ (k), which is not analytic but meromorphic in k 0 .The pole singularities also give additional contributions.
This problem is serious since the Feynman rules with energy-momentum conservation like in Eq. ( 2) are already assumed in all the approaches discussing the integration contours in the complex energy plane.This implies that those approaches have no ground in QFT.
We organize the paper as follows.In order to discuss the complex ghost problem properly in QFT, we adopt in this paper Lee's purely scalar field theory model [9], which Lee devised to mimic the essential features of the fourth-order derivative system like finite QED or QG.We will present the Lee model in Section II and its quantization as an indefinite metric QFT in the manner given by Nakanishi [28].In particular, we will explain the unfamiliar metric structure of the complex ghost field ϕ and derive two expressions for the ghost propagator; the 3d and 4d momentum expressions.In the latter, as we shall see, the integration contour of the zeroth component k 0 must take a detour around the complex poles that is much deviated from the real axis.This causes a considerable complication in the Feynman diagram computations despite the apparently covariant compact 4d-momentum expressions.
Our main concern is whether the ghost and anti-ghost are really created in the scattering processes of physical particles.In Section III, we examine the simplest process of single ghost/anti-ghost production by the scattering of two physical particles, ψ + ψ → ϕ/ϕ † .We calculate this production probability to the lowest order in three ways.First is the direct calculation of the production amplitude, given in subsection III A, which is almost trivial and simply given by the complex delta function δ c , so we give its precise definition and derive some basic properties in subsection III B. We then explicitly show that the ghost is actually produced with non-vanishing probability by the two physical particle scattering if the incident energy E is above Re √ M 2 − Im √ M 2 (lower threshold) and below Re √ M 2 + Im √ M 2 (upper threshold) with M 2 being the complex mass squared of the ghost.Although this is enough for a proof of ghost appearance, we also calculate this production probability in subsection III C by computing the imaginary part of the forward scattering amplitude of ψ + ψ → ψ + ψ with the ghost/anti-ghost intermediate line (propagator) in the s-channel.We present two methods of computation for this by using 3d-and 4dmomentum expressions for the ghost propagator in subsubsections III C 1 and III C 2. These calculations give the same result as the direct calculation in subsection III A, as a result of the optical theorem.These calculations are presented not only for checking the mutual consistency between various ways of calculation, but also for showing how the calculations of Feynman graphs should be performed in the presence of complex ghost fields and complex delta functions.
We consider a two ghost production process from a two physical particle scattering in Section IV.The direct calculation of the production probability is the simplest, but it is deferred to subsection IV B. In the first subsection IV A, we present the calculation of the forward scattering amplitude ψ+ψ → ψ+ψ at one-loop in which a pair of ghosts/anti-ghosts circulates.We present this to demonstrate how the proper calculation goes since this type of loop diagram has been discussed by many others, whose calculations often have the problem mentioned above in Eq. ( 2) already at the starting expression.The calculation using 3dmomentum expression for the ghost propagator is presented in subsubsection IV A 1, which is simpler compared with the one using 4d-momentum expression presented in IV A 2. The latter calculation is actually rather tough and lengthy and its main body has been moved to the Appendix.Instead, we add there a concise explanation for the reason why the naive Feynman rule (2) is wrong and how it should be modified, since such explanations given for many examples are buried in the lengthy calculations moved to the Appendix.
Those three ways of computation give the same result as given in subsection IV B for the production probability of a ghost pair.This result contains the complex delta function which is integrated with respect to the 3-momentum q of ghosts.Since the complex delta function indirectly depends on q through the ghost energy ω q = q 2 + M 2 as δ c (E − ω q ), the q-integration reveals a new interesting aspect of the complex delta function δ c .So, in subsection IV C, we discuss how to evaluate the q-integral and obtain the clear result that the production probability vanishes for the incident energy below the lower threshold and becomes a well-defined finite value for energy above the upper threshold.For energy E between the upper and lower thresholds, the result is divergent though the smearing of the incident energy with any finite width gives a well-defined production probability, as is usual for the distribution.
Section V is devoted to the conclusion.

II. LEE'S MODEL
In order to examine the properties of complex ghost fields in as simple a manner as possible, Lee devised in Ref. [9] a purely scalar field model which mimics the essential features of the Lee-Wick's finite QED theory in this respect.We call it Lee's model here and explain this model in a clear manner as Nakanishi presented in Ref. [28].
The system consists of three real scalar fields A, B and C in the 'photon' sector and a normal scalar field ψ in the matter sector, whose free Lagrangian is given by Here, A(x) is an analogue of the 'photon' A µ in QED that possesses a small mass δ, and B(x) is a negative metric regulator field with a large mass m accompanying the 'photon' A.
The C(x) (which mixes with B and is absent in the original Lee-Wick QED) is introduced to simulate the continuum states like lepton pairs, e + e − and µ + µ − .The mixing is represented by the last term −γ 2 BC.We are using the space-favored metric η µν = diag(−1, +1, +1, +1), so that the 'photon' A and continuum C fields, as well as the matter field ψ, are of positive metric while the regulator B is of negative metric.We also note that all the fields are real (i.e., hermitian) 6 , and that the Lagrangian ( 5) is hermitian if the mass parameters δ, m and the B-C mixing γ are real.
The BC sector of the free field Lagrangian can also be diagonalized by introducing the following complex ghost field ϕ(x): Then, we have where the mass squared M 2 for the ϕ field now takes a complex value: As shown by Nakanishi [28], the complex ghost field ϕ can be canonically quantized and is expanded as where ω q is the complex energy and the creation and annihilation operators satisfy the off-diagonal commutation relations: The Hamiltonian in the BC sector is constructed as usual from the Lagrangian (7) and is given by From the commutation relations (11) and this Hamiltonian, we see that the 1-particle ghost states |α(p yield energy eigenstates with eigenvalues ω * p and ω p respectively, and possess the following inner products and norm properties: That is, these energy eigenstates are themselves zero-norm states and have non-vanishing cross-innerproducts between two eigenstates belonging to mutually complex conjugate eigenvalues.This is a general property among energy eigenstates possessing complex eigenvalues of hermitian Hamiltonians.Indeed, for two eigenvectors |A and |B corresponding to two complex eigenvalues E A ( = E * A ) and E B ( = E * B ) of a hermitian Hamiltonian H that satisfy we have This implies that the innerproduct A|B can be non-vanishing only between two eigenvectors belonging to mutually complex conjugate eigenvalues.This further implies that any eigenvector of a complex eigenvalue necessarily has zero-norm: E|E = 0 provided that E = E * .A comment is in order here on the time evolution of complex ghost states.One may think it problematic that the ghost state |β(p) evolves as e −iωpt |β(p) with a coefficient that blows up exponentially as t → ∞ since Im(ω p ) is positive.The anti-ghost |α(p) , on the other hand, might be thought unimportant since it evolves as e −iω * p t |α(p) with coefficient decreasing exponentially as t → ∞.However, it is important to note that such exponential blowing-up or decreasing of the coefficient does not imply that the same is true for the norm of the corresponding state.Indeed, since the ghost |β(p) and its ant-ghost |α(p) have nonvanishing innerproduct only with each other, one could define their creation/annihilation operators at time t as where we recall that |B(p; t = 0) is just the state created by the Schrödinger field B = (ϕ + ϕ † )/ √ 2 at t = 0.As we will address shortly, it is important that the complex ghost appears only in this real combination (superposition) in the interaction Lagrangian.The 'photon' A(x) and the matter ψ(x) are normal positive metric fields which are expanded as usual as follows: The Lagrangian of the entire system is assumed to be of the form where it is important that the interaction Lagrangian L int depends on the ABC fields only though the combined field That is, the 'photon' A is always accompanied by the complex ghost ϕ and the anti-ghost ϕ † , so that the complex ghosts are created and annihilated always in a real superposition ϕ + ϕ † = √ 2B.Examples of the interaction Lagrangian which we employ in this paper are We consider these interactions in perturbation theory by going to the interaction picture in which the unitary time evolution operator U(t, t 0 ) is given by The fields ψ(x) and φ(x) appearing here are free fields whose explicit forms are given above.The 'photon' propagator thus always appears in the form of a φ propagator, which is given by Note, however, that despite its covariant appearance, the integration contour over q 0 is not along the real axis for the ghost propagator part ∝ 1/(q 2 +M 2 ), while it is so for the 'photon' part ∝ 1/(q 2 + δ 2 − iε) and the anti-ghost part ∝ 1/(q 2 + M * 2 ). 9n order to see this, we may calculate, in particular, the propagator of the ghost field ϕ explicitly by using the plane wave expansion (9) of ϕ(x).Recalling the definition of the T-product and using the commutation relation (11), we find 0| Tϕ(x) ϕ(y) |0 = d 3 qd 3 p (2π) 3 2ω q 2ω p θ(x 0 − y 0 )e i(qx−ωq x 0 )−i(py−ωpy 0 ) 0| α(q)β † (p) |0 + θ(y 0 − x 0 )e i(py−ωpy 0 )−i(qx−ωq x 0 ) 0| α(p)β † (q) |0 = − d 3 q (2π) 3 2ω q θ(x 0 − y 0 )e iq(x−y)−iωq (x 0 −y 0 ) + θ(y 0 − x 0 )e −iq(x−y)+iωq (x 0 −y 0 ) .(27) Note that the over-all minus sign has come from the negative norm commutation relation (11) of the ghost field.In order to rewrite this last expression of 3d integration over d 3 q into the usual 4d integration form d 3 qdq 0 by introducing q 0 variable as the q 0 -integration contour C must be the one as drawn in Fig. 1 (left) which passes below the left pole at q 0 = −ω q and above the right pole at q 0 = +ω q .With this contour we can FIG. 1. Left: The integration contour C for the ghost propagator.Right: The deformed integration contour corresponding to (78).
pick up the right pole giving the positive frequency part e −iωqt for t = x 0 −y 0 > 0, since then the q 0 integration contour can be closed by adding the half circle at infinity in the lower half plane.This is as usual, but, unlike the usual real energy case, the energy ω q for the ghost ϕ has finite positive imaginary part so that the left pole is placed much below the real axis and right pole much above the real axis, implying C becomes extremely deformed for the ghost case as drawn in Fig. 1 (left).This situation becomes opposite for the anti-ghost ϕ † , i.e., the left pole at q 0 = −ω * q is already placed above, while the right pole is placed below the real axis so that the normal integral along the real axis satisfies the required property.For the 'photon' field A, however, both poles are on the real axis so that we put the usual −iε with the mass squared to specify the integration contour properly.Thus the proper form of φ-propagator is written as where φ a denotes three component fields (A, ϕ, ϕ † ) with φ = A + (ϕ + ϕ † )/ √ 2, η a is the weight factor defined as η a = (+1, −1/2, −1/2), (−1) |a| is the norm factor (+1, −1, −1) and ω a q and q µ a are the energy and the on-shell 4-momentum defined as We use the expressions ( 29) and (30) for the φ propagator in our calculations in the following sections.The second line of (29) gives the 3d momentum form of the φ-propagator with onshell q 0 a = ω a q , which is free from the complications of the integration contour of q 0 , despite the fact that it appears non-covariant.The third line (30) gives the 4d momentum form for which the q 0 integration contour should be R, C and R when φ a represents the 'photon' A, ghost ϕ and anti-ghost ϕ † , respectively.
For general Feynman diagrams possessing many propagators, the integration contour for the energy variable q 0 for each propagator must satisfy such a requirement.When the energy conservation conditions are imposed via the vertices, those energy variables become dependent on each other, which transforms the requirements on the contours of the original energy variables into much more involved conditions of the independent energy variables.One can avoid such a complicated consideration of the integration contours if the 4-th component q 0 is not introduced and one uses only the 3d momentum variables in a similar manner to the old-fashioned perturbation theory.This method lacks manifest covariance, but the deformation of the integration contour of only the energy variables violates manifest covariance anyway even if the 4d momentum expression of propagators is used.

III. SINGLE GHOST PRODUCTION
Let us now consider the simplest process of single ghost production by the collision of two matter particles: ψ + ψ → φ.We take an interaction Lagrangian of the form and calculate the ghost production probability in two ways in this section.

A. Direct calculation of φ-production
The initial state is taken to be the two particle matter state To first order in the coupling strength f of Dyson's S-matrix, S = U(∞, −∞), with U being the operator Eq. ( 25), this initial state evolves into a single φ-particle state as where , and (• • • ) represents other states consisting of (ψ † ) 4 φ † and (ψ † ) 2 φ † .Here, since the field φ stands for the linear combination of three fields, A + (ϕ + ϕ † )/ √ 2 in Eq. ( 23), the φ † (q) |0 together with the related coefficient is understood to represent the following three terms: The hat attached to q 0 is a reminder of the fact that it changes meaning depending on the following states; i.e., q 0 = (ν q , ω q , ω * q ) for the 'photon' state a † (q) |0 , the ghost β † (q) |0 and anti-ghost α † (q) |0 , respectively.
C(x 0 ) is the regularization factor that avoids the exponential divergence difficulty for x 0 → ±∞, as was noted by Lee and Wick.For C(x 0 ) we can take, for example, a sharp cut-off that represents an interaction acting only in a finite time interval x 0 ∈ [−T, +T ]: To simplify our analytic treatment, however, we adopt the following adiabatic cut-off which corresponds to a smoother counterpart of (36) with 2T ∼ √ π/a: We multiply L int in the time evolution U operator (25) by this regularization factor and eventually take the limit a → 0 (or T → ∞) to remove the regularization.Nevertheless, it is important that the unitarity of the U operator always holds for any finite a since C(x 0 )L int (x 0 ) is hermitian at any time.Performing d 3 x integration in Eq. (34) and then d 3 q integration, we obtain where we have introduced a regularized complex delta function As the adiabatic factor is removed in the a → 0 limit, this goes to the complex delta function δ c (z) with complex argument z, which in turn reduces to the usual Dirac's delta δ(x) for real argument z = x, and was first introduced by Nakanishi [26,27] long ago in connection with the construction of Hamiltonian eigenstates of unstable particles with complex eigenvalue.We shall show in the next subsection that δ c (z) gives a well-defined distribution and write henceforth only the expressions for the limit a = 0 unless the explicit regularization is necessary.
2. Left: Graphical presentation of the norm (41).Right: The forward scattering amplitude of two physical particles.
Noting Eq. ( 35), using the commutation relation (11) as well as [ a(p) , a † (q) ] = δ 3 (p−q) and writing p µ 1 + p µ 2 = P µ , we can calculate the norm of the produced φ-state, (S − 1) (1) |I φ in Eq. (38), as where we have distinguished the momentum P ′ of the bra-state P ′ | from that P of the ket-state |P for clarity although P ′ = P when computing the norm |P 2 = P |P .
Fig. 2 (left) is a graphical presentation of the norm (41).We have also used the identity which is shown to hold for the complex delta function δ c (P 0 − q 0 ) with real P 0 and complex q 0 , shortly below.
In the norm expression (41), the factor (2π) 3 δ 3 (P ′ − P )(2π)δ(P 0 − P 0 ) = (2π) 4 δ 4 (0) is the total space-time volume so that it is divided out for the production probability per unit space-time volume, as usual.The first term is the norm of the 'photon' particle A and the second and third terms are that of the ghost and anti-ghost, respectively, with negative sign.The first term yields the production probability P A of 'photon' A per unit space-time volume as This implies that the 'photon' production occurs only exactly at the resonating energy P 0 = ν P = √ P 2 + δ 2 while the production probability is zero below and beyond the threshold P 0 = ν P .This term merely reproduces the usual result, confirming the validity of this computation.
The second and third terms, therefore, give the production probability P ϕ of our ghost (superposition) ϕ + ϕ † as We note that this probability is real and negative as shown more explicitly in the next subsection, thus implying a violation of physical unitarity: That is, if it is non-vanishing, the probability of other final physical states (consisting of physical particles alone) exceeds one.
Before ending this subsection, we add an important remark here: The opposite process to the single ghost production is the "decay" of the ghost into two physical particles, which of course occurs due to time reversal invariance.One may then naturally wonder, are the complex ghost particles stable?The answer is yes (in contrast to the assumption of [17,18]).First of all, the ghost state created in the superposition ϕ+ϕ † has a negative norm.Since the Dyson's S-matrix of the present system is unitary, the negative norm, say −1, of the initial ghost state must be conserved.So, whatever final states are produced from the initial ghost state, the norms of all those final states sum up to the value −1 of the initial ghost state's norm.To realize this negative value, however, ghost particles must be contained among the final states.This implies that the ghosts can never disappear by totally 'decaying out' into lower mass physical particles.

B. Complex delta function
To see that P ϕ is actually non-vanishing, let us now analyze the property of the complex delta function δ c (P 0 − ω P ) in more detail.From the definition (40), we have an expression where the ghost energy ω P is the complex quantity In the center of mass (CM) frame, P = 0, this ω P becomes The ∆ a function then becomes (denoting P 0 simply by E) This function, as a function of (real) energy E, is clearly non-vanishing everywhere for finite a, so that the ghost production probability P ϕ in Eq. ( 44) is non-vanishing for any initial energy E = P 0 on any finite time interval T ∼ a −1 .This may not sound surprising since the energy is not conserved for finite time interval in any case.Now consider the infinite time limit a → 0. Note that the parabolic function of E in the exponent in Eq. ( 48) Inside this interval, on the other hand, the function ∆ a (E − ω) is non-vanishing but diverges while rapidly oscillating as a → 0, implying that the limiting function δ c (E − ω) is not a function, but rather a distribution.To see that it gives a well-defined distribution possessing non-vanishing support in the interval |E − Re ω| < Im ω, let us compute an integral of it multiplied by a test function f (E) which is analytic in a necessary domain: We can deform the contour of this E integral along the real axis R to the following: where the deformed contour is shown in Fig. 3.The integrals along the finite vertical segments kept finite.The integral can thus be evaluated by the contour integral along R(ω), which denotes the horizontal contour parallel to the real axis R and passing the point z = ω i.e., z = x + iIm ω (with real x ∈ [−∞, ∞]).This deformation is allowed when the test function f (E) is analytic in the complex plane inside the rectangular region surrounded by R + C 1 + R(ω) + C 2 (see Fig. 3).The integral along R(ω) is evaluated by making the change of variable which proves that the distribution δ c (E − ω) works as if it is the usual Dirac's delta function despite the fact that E is real and ω is complex.This also proves which proves Eq. (42) when f (E) = δ c (E − ω), as promised.
We should note that the remarkable property (54) for the complex delta function holds if and only if the test function f (E) has no singularity inside the rectangular region stated above [28].In practice however, f (E) is usually given by the products of meromorphic Feynman propagators which often have pole singularities in the rectangular regions in question, in which case one has to pick up that contribution also.As we will see soon, this is actually the point touching the core of the problem.
In the actual scattering experiment, the total energy E of the initial particles, P 0 = p 0 1 +p 0 2 , must necessarily have a certain uncertainty around P 0 , which may be described, for instance, by a Gaussian distribution with standard deviation σ of the form For such initial particles, the ghost production probability (44) should be averaged with this distribution.The actual production rate of complex ghost per unit space-time volume is then found to be which is finite as far as σ is finite.

C. Imaginary part calculation of the forward scattering amplitude
Since the interaction Hamiltonian is hermitian, the total S-matrix unitarity follows trivially and leads to the optical theorem for the T -matrix: where T |I F denotes all independent states |F contained in the final scattered state T |I .
In perturbation theory, this optical theorem holds for each set of the same order terms in the coupling constant f .Additionally, the imaginary part of the particular Feynman diagram on the LHS equals the RHS with the sum over F restricted to a certain subset.
Let us now calculate the forward scattering amplitude of two physical particles with the initial state |I(p 1 , p 2 ) in Eq. (33).At second order in f of the interaction Lagrangian (32), there are only three diagrams which have φ-propagator in s-, t-and u-channels, and we calculate only the s-channel amplitude since it has an imaginary part that reproduces the φ-production probability in the previous subsection III A.Here we note that the result Eq. (41) which we have found as the norm of the produced φ-particle state (34) in the subsection III A can be rewritten in the form of the RHS of optical theorem (58) as (2π) 4 q δ 4 (P F − P I ) T |I(p 1 , p 2 ) φ(q)-prod.
We should find the same result by calculating the imaginary part of the forward scattering amplitude on the LHS of the optical theorem (58).The s-channel tree level amplitude at O(f 2 ) corresponding to the Feynman diagram Fig. 2 (right) is now also easily computed: where P µ := p µ 1 + p µ 2 , P ′µ := p ′µ 1 + p ′µ 2 and C a (x) is the adiabatic Gaussian cutoff (37).We will compute this S-matrix element (60) in two ways; first using the 3d-momentum form (29) of the φ propagator, and second using the covariant 4d-momentum form (30).We will explicitly see why a careful treatment of a regularization (cutoff) (37) is indispensable for obtaining a coincident final result in both calculations.Once the equivalence of these methods is established, we can use either form of the propagator freely.

Use of 3d-momentum form of the φ propagator
We first present the calculation using the 3d-momentum form (29) of the φ propagator.Now, we insert the 3d-momentum form (29) into the propagator 0|Tφ(x)φ(y)|0 in Eq. (60).Then, making a change of time variables (x 0 , y 0 ) into CM and relative time (T, t) and writing explicitly the adiabatic (Gaussian) cutoff factor (37) at both interaction points x 0 and y 0 , we have, after performing the d 3 xd 3 y integrations and using δ 3 (P ′ − q)δ 3 (q − P ) = δ 3 (P ′ − P )δ 3 (q − P ), LHS of (60) × dT e −2a 2 T 2 e i(P ′0 −P 0 )T dt e −a 2 t 2 /2 θ(t)e i(P 0 −ω a q )t + θ(−t)e i(P 0 +ω a with P 0 := (P ′0 + P 0 )/2.Note that the integrations of the CM time T and relative time t are totally separated including their adiabatic cutoff factors.The dT integration yields the factor dT e −2a 2 T 2 e i(P which together with the above (2π) 3 δ 3 (P ′ − P ) completes the total energy-momentum conservation factor i(2π) 4 δ 4 (P ′ − P ), which should be divided out when giving a T -matrix element.
Next, we consider the dt integration part (denoting P 0 by E as in (48)) For the anti-ghost case first, ω a q = ω * q has negative imaginary part Im ω * q < 0. We note that the integral converges even in the limit a → 0 since e i(E∓ω * q )t decreases exponentially as t → ±∞, meaning that the a → 0 limit is easily calculated by setting a = 0 from the beginning: (66) In the ghost case, ω a q = ω q has positive imaginary part Im ω q > 0 oppositely to the antighost case.So the a → 0 limit is divergent since e i(E∓ωq)t exponentially blows up as t → ±∞, respectively.But this also implies that the functions e i(E∓ωq)t decrease exponentially in the opposite directions t → ∓∞.If we rewrite θ(±t) as 1 − θ(∓t) in Eq. ( 65), then those terms with a −θ(∓t) part give convergent quantities in the a → 0 limit as in the anti-ghost case, while the terms with 1 give the regularized complex delta function; that is, we can evaluate Eq. (65) in the ghost case as lim a→0 dt e −a 2 t 2 /2 θ(t)e i(E−ωq)t + θ(−t)e i(E+ωq)t = lim a→0 dt e −a 2 t 2 /2 e i(E−ωq)t + e i(E+ωq)t − dt θ(−t)e i(E−ωq)t + θ(+t)e i(E+ωq)t = lim Finally, we consider the 'photon' which has a real energy ω a q = ν q .In this case, we can follow the usual −iε trick to shift the mass square which results in the energy shift ν q → ν q − iε.Then, since Im (ν q − iε) < 0, this reduces to the above anti-ghost case and the result (66) immediately gives dt θ(t)e i(E−νq+iε)t + θ(−t)e i(E+νq−iε)t = i If we recall the well-known formula 1/(x ∓ iε) = P(1/x) ± iπδ(x) (with P denoting principal value), this result can be rewritten into the form Eq. (68 It is interesting to note that this can further be rewritten into which is the same result as the one which we would have obtained if we had used +iε trick to replace ν q → ν q + iε and applied the formula (67) for the ghost case.Since (68) equals (70), both tricks replacing ν q by ν q ± iε lead to the same result.Now, let us go back to the Eq.(63).We apply these formulas (68), ( 67) and (66) to the dt integration part (65) for φ a = 'photon' A, ghost ϕ and anti-ghost ϕ † , respectively, and divide out the total energy-momentum conservation factor i(2π) 4 δ 4 (P ′ − P ), setting P = P ′ in the Eq. ( 63), to find the forward scattering transition amplitude where we have used ) in the presence of δ 3 (P − q).We thus finally obtain the imaginary part of the s-channel forward transition amplitude, where we have omitted the terms δ(P 0 + ν P ) and δ c (P 0 + ω P ) since they vanish for P 0 > 0. 10This imaginary part exactly coincides with the RHS result (41) divided by (2π) 4 δ 4 (0), that gives the φ-production rate per unit space-time volume directly calculated in the previous subsection III A, thus confirming the optical theorem (58).

Use of 4d-momentum form of the φ propagator
Now we present the calculation of the S-matrix element (60) using the 4d-momentum form (30) of the propagator: where we recall that the q 0 -integration contour for the ghost ϕ is the much deformed C drawn in Fig. 1 (left), while those for 'photon' A and anti-ghost ϕ † are along the real axis R.Then, performing the d 4 x and d 4 y integrations yields11 and after dividing out (i times) the total momentum conservation factor, i(2π) 4 δ 4 (P ′ − P ), and setting P ′ = P , we find the following expression for the T -matrix element for the forward scattering s-channel diagram: Suppose that we perform the d 4 q integrations in this Eq.(75) to simply put q µ equal to P µ by naively using the complex delta function δ c (q − P ), as if it were just like the Dirac delta function.We would then immediately obtain 12   I| T |I s-ch.
which is exactly the expression that one would write down if the usual Feynman rule in momentum space is applied directly to the diagram Fig. 2 (right), in which the naive energymomentum conservation at each vertex is tacitly assumed in advance.This represents a very elementary, but very serious mistake which many have made, including Lee and Wick, and even many of those who were critical of their theory for other reasons.If Eq. ( 76) were correct, the ghost contributions, the second and third terms, add up to a real quantity and thus do not contribute to the imaginary part of the transition amplitude.Only the first 'photon' part gives the following imaginary part (for the s-channel P 0 > 0): where the term proportional to δ(P 0 + ν P ) = 0 is suppressed.Essentially by this type of calculation and reasoning, Lee and Wick concluded that the complex ghosts are not produced from physical particle scattering.The correct computation of the RHS of Eq. ( 75) should be performed as follows, noting that the naive use of the complex delta function δ c (q 0 − P 0 ) is allowed when the integration variable q 0 remains real all the way along the contour since then δ c (q 0 − P 0 ) is reduced to the usual Dirac delta function δ(q 0 − P 0 ) with real argument q 0 − P 0 .Since this is already the case for the first term ('photon' part) and the third term (anti-ghost part), the problem lies in the second term (ghost part), for which the integration contour C is much deformed one as depicted in Fig. 1 (left) on the complex q 0 plane.To avoid the moving imaginary part Im q 0 on C, we can deform the contour C according to as shown in Fig. 1 (right), where C(ω) denotes an infinitesimal circle rotating anti-clockwise around the pole at q 0 = ω.Then, the first integral part R dq 0 is again along the real axis, for which we can use the naive substitution rule q 0 → P 0 also for the second term (ghost part).This contribution together with those from the first and third terms in (75) reproduces Lee and Wick's amplitude (76), which results from applying the usual Feynman rule with energy-momentum conservation assumed at each vertex.
Here, however, we have two extra contributions from the two circles C(±ω q ) which can be evaluated by the pole residues according to Cauchy's theorem: dq 0 1 −q 0 2 + ω 2 q δ c (q 0 − P 0 ) 12 The sum of the ghost and anti-ghost propagators in (76) is the real propagator for the fakeon [15,16]. 13It is important to note that the complex delta function δ c (q 0 − P 0 ) is the a → 0 limit of the analytic (Gaussian) function ∆ a (q 0 − P 0 ) which has no singularity in the relevant domain for the deformation performed in Eq. (78).
Adding this to the Lee-Wick's amplitude (76), exactly reproduces the previous result (71) obtained by using the 3d-momentum propagators.
At this stage we emphasize that in the 3d form of the propagators (27), the time propagation of the positive and negative energy components are clearly distinguished thanks to θ(±t) function, where, in the case of the ghost, we mean a positive or negative energy with respect to its real part.As a result, there is no violation of causality even in the presence of the ghost.However, the amplitude (76) will violate causality [11,17,18] (see also [8]), because the integration contour of the 4d form of the ghost and anti-ghost propagators (which inherits the causal information from the 3d form of the propagators) is completely ignored.

IV. GHOST PAIR PRODUCTION
The one-loop amplitude for ψ+ψ → ψ+ψ shown in Fig. 4 (right) has been representatively considered to demonstrate the Lee-Wick prescription [8][9][10] of how to integrate the loop momenta.According to their prescription, the imaginary part of the amplitude should vanish for the parts where ghost and/or anti-ghost are propagating in the loop.The optical theorem would then imply that the production of two ghost particles by a collision of two physical particles, i.e. ψ + ψ → B + B (B = (ϕ + ϕ † )/ √ 2) is forbidden.Here in this section we will examine their calculation in our field theory framework.To this end, we consider the interaction Lagrangian A. Imaginary part calculation of the forward scattering amplitude 1. Use of 3d-momentum form for the φ propagator Since we are interested in ghost production, we focus only on the s-channel contribution to the one-loop amplitude for the scattering ψ+ψ → ψ+ψ.With the interaction Lagrangian L int of (80), the calculation proceeds quite in parallel to that in the previous section.As in Eq. (60), we have at second order in the coupling f with the initial 2ψ state |I(p 1 , p 2 ) defined in (33) and , where C a (x 0 ) is the adiabatic cutoff introduced in Eq. (37) and Σ(x − y) is given by Here, φ a = (A, ϕ, ϕ † ) are the three component fields in φ = A + (ϕ + ϕ † )/ √ 2 with the weight factor η a = (1, −1/2, −1/2) and norm sign (−1) |a| = (+1, −1, −1).Inserting the 3d momentum expression in Eq. (29) for the propagator 0| T φ a (x)φ a (y) |0 , we find a compact form for Σ ab : where ω a q and q µ a are the energy and the on-shell 4-momentum of φ a defined before in Eq. ( 31).Noting the similarity between Eqs. ( 60) and (81), and recalling the procedure to obtain Eq. ( 63), we perform the d 3 xd 3 y integrations in (81), yielding three dimensional delta functions (2π) 3 δ 3 (q + q ′ − P ) and (2π) 3 δ 3 (q + q ′ − P ′ ).This then allows us to carry out the integration over q ′ to write q ′ = −q + P = −q + P ′ . 14Then, as we have done in (61), we use the CM and relative time (T = (x 0 + y 0 )/2, t = x 0 − y 0 ) to perform the time integration and obtain = −f 2 (2π) 3 δ 3 (P − P ′ ) dT e i(P 0 −P ′0 )T e −2a 2 T 2 14 The second term in (83), i.e., the one proportional to θ(−t), actually gives q ′ = −q − P = −q − P ′ , but we make the change q → −q which implies −q − P → q − P .We then use the fact that ν q = ν −q and dt e −a 2 t 2 /2 θ(t) e i( P 0 −q 0 ab )t + θ(−t) e i( P 0 +q 0 ab )t , where P 0 = (P 0 + P ′0 )/2 = (p 0 1 + p 0 2 + p ′0 1 + p ′0 2 )/2 and q0 ab = ω a q + ω b P −q , and the T integration gives (2π)δ(P 0 − P ′0 ) in the a → 0 limit.When applying the formulas to carry out the t integration with either (67), (66), or (68) and θ(±t), we must classify the cases when the imaginary part of q0 ab = ω a q + ω b P −q becomes positive or negative, which is a bit complicated task.In order to avoid such an inessential complication for the present issue of ghost production, we discuss the problem in the CM frame and set P = 0 henceforth.Then, the equality ω a P −q = ω a q holds so that we have Eq.( 68) Im q0 ab > 0 : (A, ϕ), (ϕ, ϕ) Eq.( 67) Im q0 ab < 0 : (A, ϕ † ), (ϕ † , ϕ † ) Eq.( 66) . (85) Applying the indicated formula for the dt integrations in (81), dividing out the factor i(2π) 4 δ 4 (P ′ − P ) and setting P = P ′ , we obtain the forward scattering T amplitude where we have omitted the terms δ c (P 0 + 2ω q ) and δ c (P 0 + ν q + ω q ) which vanish for the s-channel P 0 > 0 because of (51), and used P 0 = (P 0 + P ′ 0 )/2 = P 0 .Note that we applied Eq. ( 68) for the (ϕ, ϕ † ) loop because Im (ω q + ω * q ) = 0.This result (86) for the one-loop T amplitude consists of fraction and complex delta function pieces.All the fraction parts other than that from the (ϕ, ϕ † )-loop (containing (ω q + ω * q ) 2 ) are identical with the result of Lee [9].However, the complex delta functions and the fraction part from the (ϕ, ϕ † )-loop with −iε, which contribute to the imaginary part corresponding to the ghost production rate, are all absent in Lee's computation [9].
From our result (86), the correct imaginary part of the forward scattering T amplitude is found to be 2Im Again, we have omitted all the terms of the form δ c (P 0 + ω a q + ω b q ) which vanish when P 0 > 0. Needless to say that only the first term ('photon-photon' term) is present in Lee's computation.
These delta functions, other than the last one, are all the usual Dirac functions.Using these usual delta functions we can carry out the q ′ integration, divide out the factor i(2π) 4 δ 4 (P ′ − P ) and set P = P ′ to obtain the forward scattering T amplitude where q ′ = P − q is understood, and ω a q = (ν q , ω q , ω * q ).We should however emphasize that the q 0 -or q ′0 -integration using the complex delta function δ c (P 0 − q 0 − q ′0 ) is very non-trivial because the argument P 0 − q 0 − q ′0 is complex when φ a or φ b is the ghost ϕ.In such cases, we have to make a shift of the integration contour to make q 0 + q ′0 real so as to reduce the complex delta function δ c to the usual Dirac delta function δ, thus allowing it to pick up the pole singularities of the propagators.
We may then perform this task for all the 3 × 3 = 9 cases of combination of two field propagators in Eq. ( 90) systematically.This task is straightforward but becomes slightly tedious and a bit lengthy, so we move the calculation to the Appendix.There we actually obtain exactly the same result as the above Eq.(86) obtained by using the 3d form of propagators.
Here, however, we should concisely explain the reason why the usual Feynman rule expression (2) with energy-momentum conservation used in advance is wrong and does not follow from the correct expression (3), as announced in the Introduction.The proofs for this are given for many concrete cases in the Appendix, however, since the essential point is buried in the lengthy calculation, it would be better to recapitulate it here.
We start with the correct expression (3) in the Introduction, which reads, after performing the 3d k integration by using the usual Dirac delta function δ 3 (k + q − p), where ω q = q 2 + M 2 and E k = k 2 + µ 2 are the energies of the complex ghost ϕ and physical particle ψ, respectively.We have written the k 0 integration contour as R (real axis) explicitly while putting the usual −iε to indicate how to avoid the pole on the real axis.This expression corresponds to the case for φ a = ϕ (ghost) and φ b = ψ (physical particle in place of the photon A) in Eq. ( 90).Now the problem is the fact that q 0 is complex since the q 0 -integration contour C is the much deformed contour in the complex q 0 plane as shown in Fig. 1.So, when using the complex delta function δ c (k 0 + q 0 − p 0 ) for k 0 integration, we need to make the argument k 0 + q 0 real by shifting the integration variable k 0 in order to reduce the δ c (k 0 + q 0 − p 0 ) to the usual Dirac delta function.To do so, it is better to make the imaginary part Im q 0 not run.So we change the q 0 -integration contour C in Fig. 1 (left) to the contour Fig. 1 (right) as we did in Eq. (78) for the single ghost production case.Then, the q 0 -integration part in Eq. ( 91) can be rewritten into The last ± terms are the contributions from C(±ωq) dq 0 around the poles at q 0 = ±ω q .Note that the complex delta function δ c (k 0 + q 0 − p 0 ) has become the usual Dirac delta function δ(k 0 + q 0 − p 0 ) in the first term on the RHS since q 0 remains real on the contour R. Inserting this and performing the trivial k 0 integration using this Dirac delta function in the first term, the Eq.(91) becomes Eq.(91 The k 0 -integration for the second two (±) terms is non-trivial.To make the arguments k 0 ± ω q − p 0 of δ c real there, we shift the k 0 -integration contour from the real axis R to the contour R(∓ω q ), respectively.Here R(z 0 ) with a complex z 0 generally indicates the horizontal contour parallel to the real axis and passing the point z 0 .However, the point is that this shift of the k 0 -integration contour from R to R(∓ω q ) does not change the integral, only if there are no singularities in the k 0 integrand in the k 0 -domain between the two contours R and R(∓ω q ) (provided that the contributions from the vertical passes at ±∞ between them vanish as is the case here).In this rectangular domain, however, the present integrand 1/(E 2 p−q − k 0 2 − iε) actually has a pole at k 0 = ±(E p−q − iε), respectively, so that this shift also picks up the contribution from the poles and we have with C ±(E p−q − iε) denoting the infinitesimal circle surrounding the pole k 0 = ±(E p−q − iε) anti-clockwise.Note that k 0 ± ω q is real on the contour R(∓ω q ) for the first term so that δ c (k 0 ± ω q − p 0 ) reduces to the Dirac delta function and hence the k 0 integration becomes trivial.Noting also that the integral for the second term is evaluated by the Cauchy theorem, we find Substitution of this into the second term in Eq. (93) gives Eq. (91 We immediately notice that the first plus second terms in fact give the original usual Feynman rule expression (2) with normal energy-momentum conservation used: Note that the integration contour is C.So we find that this naive Feynman rule misses the extra contributions expressed by the third term in Eq. ( 96).
If combined with the contribution from the ϕ † -ψ loop, the naive Feynman rule gives only a real amplitude as Lee [9] showed even if the deformed contour C is used correctly.The imaginary parts solely come from the extra third term in Eq. (96), which originates from the complex delta function.

B. Direct calculation of two φ production
Because of the optical theorem (58), the imaginary part of (81) should be equal to the norm of the intermediate state that also enters the unitarity sum: where, to the first order in the coupling f of L int in Eq. ( 80), and C a (x 0 ) is the adiabatic cutoff given in (37).In this expression (99) we have omitted the other states corresponding to the disconnected diagrams which are irrelevant to the present discussion.Further, ω a q , |η a | and φ † a (q) with index a running over three component fields φ a = (A, ϕ, ϕ † ) are given by respectively.The integration over x gives a three-dimensional delta function (2π) 3 δ 3 (P − q − q ′ ) with P = p 1 + p 2 , which can be used to perform the integration over q ′ .Then, performing also the x 0 integration by using the definition of the complex delta function (40), we obtain in the a → 0 limit We note that is just the superposition of one-particle states created by the field φ(x) = A + (ϕ + ϕ † )/ √ 2 at time x 0 = 0 as introduced in Eq. ( 35), whose norm has been already calculated when obtaining Eq. (41).In place of the norm, the inner-product is also computed in the same way as Since the state appearing here (101) is just a two-particle (tensor product) state of (102), the norm of (101) can immediately be found from Eq. (103) to read Dividing out the factor (2π) 4 δ 4 (0) and writing the sum over a, b explicitly, we find the following expression for the norm of the produced state T |I in the CM frame (P = 0): We see that this result giving the RHS of the optical theorem (58), exactly coincides with the imaginary part of the forward scattering amplitude given in Eq. ( 87) which we calculated in the previous subsection.The first term in the first line gives the production probability of two 'photon' A-A which is of course positive as usual.The second term in the first line and the two terms in the third line give the production probability of two ghosts via ϕ-ϕ † , ϕ-ϕ and ϕ † -ϕ † , respectively, and are positive due to the two negative norm particles.The two terms in the second line give a negative probability for the production of a 'photon' and a ghost, A-ϕ and A-ϕ † .These probabilities are all, except for the normal two 'photon' case, proportional to the complex delta function.We have already shown in Section III that the complex delta function is well-defined and non-vanishing, so these probabilities for two ghost particles are clearly non-vanishing and violate the physical particles' unitarity.Since these probabilities contain the 3d momentum integration d 3 q of the complex delta function, this shows a new interesting aspect of this complex delta function, so that we may now discuss more explicitly the production probability of two ghosts, ϕ-ϕ † , ϕ-ϕ and ϕ † -ϕ † .We thus have explicitly confirmed that the optical theorem (58) is satisfied for our scattering process in the lowest non-trivial order.This theorem is a trivial identity which directly follows from the unitarity of the Dyson's S-matrix as a result of the hermiticity of our interaction Lagrangian.So although the confirmation itself of the theorem does not have so important meaning, it gives a useful consistency check for the validity of our computations.
C. Explicit evaluation of the ghost production probability

Two ghosts production
Let us now investigate more explicitly the production probability for two ghosts, ϕ-ϕ † , ϕ-ϕ and ϕ † -ϕ † , which were given in the second term of the first line and the two terms in the third line in Eq. (104), respectively.
First consider the ghost-anti-ghost ϕ-ϕ † pair production, whose probability is given by This pair production is very special since the total energy ω q + ω * q of the complex ghost pair can be exceptionally real in the CM frame with P = 0.The delta function δ(ω q + ω * q − P 0 ) here is the Dirac delta function implying the usual energy conservation which determines the magnitude of momentum |q| ≡ q of produced ghost as where the total energy is P 0 = 2E.It is interesting to note that the imaginary part γ 2 of the complex ghost mass squared M 2 = m 2 + iγ 2 effectively contributes γ 4 /4E 2 to the real mass squared.Anyway, the integration over q in Eq. ( 105) can be done by the usual formula for the Dirac delta function as Next consider the 'ghost-ghost' and 'anti-ghost-anti-ghost' production whose probability is given as follows by again writing P 0 = 2E and using the property δ c (2z) = (1/2)δ c (z): In subsection III B, we have considered some properties of this complex delta function δ c (E − ω q ) as a function of E with a fixed complex ω q .Here, we have to integrate this function over the variable q = |q| and need its property as a function of q, though we have no convenient formula for the complex delta function as we do for Dirac's delta, i.e., This is the formula for the real function f (x) of real variable x, which we have just used above.We shall show a similar formula holds for some cases of the present complex delta function δ c (E −ω q ).In order to treat this problem properly, we use the adiabatic (Gaussian) regularization form of the complex delta function, so we discuss where we recall that Here, the real part of I(E) multiplying the const., i.e., (f 2 /16π) × 2Re [I(E)], gives the above 'ghost-ghost' and 'anti-ghost-anti-ghost' production probability (108).
As noted in Sect.3, the real part G R (q) of the exponent factor G(q) = (E−ω q ) 2 determines the magnitude of the regularized complex delta function ∆ a (E − ω q ).In particular, in the a → 0 limit, its sign is crucial since the form ∝ exp − G R (q)/4a 2 implies that when G R (q) > 0 the limit vanishes while it diverges (and rapidly oscillates) if G R (q) < 0. That is, the complex delta function (distribution) δ c (E − ω q ) defined as lim a→0 ∆ a (E − ω q ), when viewed as the function of E for each fixed q, has non-vanishing support only in the negative G R (q) region.Since G R (q) is quadratic in E, the negative G R (q) region is simply given by E − (q) < E < E + (q), E ± (q) = Re ω q ± Im ω q . (112) We can plot the positive/negative G R (q) region, or the lower boundary curve E = E − (q) and also the upper boundary curve E = E + (q) as a function of q in the real (q, E) plane.See Fig. 5 in which these functions are plotted for the choice of parameter γ 2 /m 2 = 0.5.As it FIG. 5.The unshaded region is the negative G R (q) region, i.e., where G R (q) < 0, and thus gives the support of the complex delta function δ c (E − ω q ).The upper and lower boundary curves of the negative G R (q) region are given by E = Re ω q ± Im ω q ≡ E ± (q).The curve in the middle of the region is E = Re ω q drawn for reference.This figure is drawn with parameter γ 2 /m 2 = 0. is easily shown, Re ω q (Im ω q ) is a monotonically increasing (decreasing) function of q(> 0) and both curves E = E ± (q) are monotonically increasing.We have called the boundary energies E − (0) and E + (0) at q = 0 the lower threshold E Lth and the upper threshold E Uth , respectively, in Sect.3.They are the lowest points of the lower and upper boundary curves E = E ∓ (q), respectively.
"contour R + (q 0 )" for a reason which becomes clear shortly.On the complex z = q 2 plane, R + (q 0 ) is a straight line parallel to the real axis (though on the complex q plane, it looks like a part of a square root curve of course).So, for the evaluation of the original integral in Eq. ( 110), we should make the following contour deformation on the complex z-plane: We can forget the contribution from the last vertical line C 2 at Re z = ∞ which vanishes as usual.And, as we will do shortly, the contribution from the horizontal contour segment R + (q 0 ) can be very easily evaluated since E − ω q is real on the whole contour R + (q 0 ) and hence, the complex delta function δ c (E − ω q ) is reduced to the usual Dirac delta function δ(E − ω q ) of real variable x.
Before doing this task, let us examine the contribution from the first vertical line C 1 in which a critical difference appears between the cases with energy in Region II and Region III.

Region II: E
As we see in Fig. 5, the q = 0+ point in Region II is in the negative G R (q) (unshaded) region.This must also be the case even when q is extended to be complex; around the origin q = 0 in the complex q plane, the real part of G(q), G R (q), is negative.This means that the complex delta function δ c (E − ω q ) = lim a→0 ∆ a (E − ω q ) becomes divergent (and rapidly oscillating) in the a → 0 limit in a finite neighborhood of q = 0. 15So, the a → 0 limit of the integral (110) is already divergent in the contribution from the contour segment C 1 alone.Therefore, there is no particular reason to adopt the deformed contour (115).We use the original integration contour R + .
Although it is divergent, the integral I(E) has a well-defined meaning as a distribution.We average the integral (110) by using the Gaussian smearing function (56) around an energy E in Region II with standard deviation σ leading to This now gives a well-defined finite value thanks to the finite width of σ although there is some interval 0 ≤ q < q + where Re [(ω q − Ē) 2 ] < 0 for the present energy Ē in Region II.
Region III: E > E Uth Contrary to the above energy Region II, the Fig. 5 shows that the q = 0 point in Region III is in the positive G R (q) (shaded) region.Again, by analyticity, this must be the case around q = 0 also in the complex q-plane or complex q 2 = z plane.Actually, for energy E in this Region III, we can easily show that the whole vertical line C 1 is in the positive G R (q) region as follows.
On the vertical line C 1 , the complex variable z is parametrized as z = 0 − iy by the real parameter y ∈ [0, γ 2 ] and ω q looks like ω q = m 2 + i(γ 2 − y).As we go down from the origin z = 0 to the point z = −iγ 2 touching R + (q 0 ), the real parameter y changes from 0 to γ 2 and the real part G R (q) of G(q) = (E − ω q ) 2 monotonically increases: Since it is positive at the starting point z = 0 for the energy E in the Region III, the G R (q) keeps positive all the way along C 1 .
The positivity of G R (q) on C 1 means that the contribution from the segment C 1 to the integral (110) vanishes in the a → 0 limit.Thus the contribution can come only from the horizontal line R + (q 0 ).Now, ω q is real on R + (q 0 ) as stated above, so G(q) = (E − ω q ) 2 is positive aside from the possible isolated zeros.The zeros, say q = q 0 , of G(q) are also easily found: That is, aside from the sign ±, the zero q 0 is uniquely given by We denote the zero with Re q 0 > 0 as q 0 and the other as −q 0 .The one on R + (q 0 ) is q 0 , so the name R + (q 0 ) means that it passes through the point q 0 .We now know that the complex delta function δ c (E − ω q ) becomes the usual Dirac delta function if written in terms of the real variable x = z+iγ 2 = q 2 +iγ 2 on R + (q 0 ) and it has only one zero x 0 = q 2 0 +iγ 2 = E 2 −m 2 .We can apply the formula (109) to δ c (E − ω q ) = δ f (x) with f (x) = E − ω q to evaluate the contribution of this contour segment R + (q 0 ) to the integral (110).
It is truly remarkable that such a non-trivial integral I(E) of the complicate function (q 2 /ω 2 q )δ c (E − ω q ) of q including complex delta function, which is actually a distribution possessing non-localized support, can be evaluated analytically by computing a single complex root q 0 of G(q) = 0.However, to verify the equivalence of the original integration over R + with that over the deformed one (115), we must confirm that the function q 2 /ω 2 q multiplied by the complex delta function δ c (E − ω q ) in Eq. ( 110) has no pole singularities in the rectangular domain surrounded by R + and C 1 + R + (q 0 ) + C 2 in the z plane.The only singularity is the pole at ω 2 q = 0, i.e., at z = −m 2 − iγ 2 which is just on the negative side line R − (q 0 ) of our integration contour segment R + (q 0 ), thus indicating clearly that the pole is outside the rectangular domain.
In Fig. 6 we compare the exact result (120) with a numerical calculation, where we have used: m = 1 , γ = 1/ √ 5 with a finite cutoff a = 1/14.The blue line presents the numerical result of the real part of (110), and the red line shows the real part of the exact result (120), which is applicable starting at E = E Uth = 1.272.We see an excellent agreement of the numerical result with the exact result already at finite a = 1/14.

Photon and ghost production
The production probability of photon and ghost, A-ϕ and A-ϕ † , is given by the second line of (104): We therefore consider the integral where √ π e − G(q)/4a 2 , G(q) = (E − (ν q + ω q )/2) 2 , ν q = q 2 + δ 2 , ω q = q 2 + m 2 + iγ 2 , GR (q) = Re G(q) = E − (ν q + Re ω q + Im ω q )/2)(E − (ν q + Re ω q − Im ω q )/2 , Comparing the expressions in (123) with those of (111) we can easily convince ourselves that the same classification of the energy regions exists here as in the case of two ghosts production: The lower and upper threshold energies are given by E Lth = (ν 0 + Re ω 0 − Im ω 0 )/2 and E Uth = (ν 0 + Re ω 0 + Im ω 0 )/2, respectively.In Region I (E ≤ E Lth ) the real part GR (q) is positive for ∀q > 0, so that the integral I ph (E) given in Eq. ( 122) vanishes in the a → 0 limit.In Region II the integral diverges, which means that we have to average it with the Gaussian smearing function (56) to obtain a finite result.The situation in Region III is slightly different: There no longer exits such a simple variable like z = q 2 in Eq. (114) that the line, realizing real ν q + ω q in the complex plane of that variable, becomes a straight line parallel to the real axis.Fortunately, it is not a necessary condition for the integral (122) to be analytically evaluated by computing a single complex root q 0 : It is sufficient for this if there exists a curve (contour) on the complex q plane, which starts from the origin (i.e., q = 0), goes through the zero q 0 of (ν q + ω q )/2 − E and approaches +∞, such that the real part of G(q) is all the way positive, except at q = q 0 .An example is shown in Fig. 7, where we have used representative values E/m = 0.9 , γ/m = 0.5 , δ/m = 0.5 with m = 2.As far as such a contour exists, the calculation of the integral reduces to the evaluation of the contribution only in the infinitesimal neighbourhood of q 0 .We thus expand the argument of the complex delta function around q = q 0 : (ν q + ω q )/2 − E = (q − q 0 )A + O((q − q 0 ) 2 ) , where (ν q + ω q )/2 − E q=q 0 = 0 and A = d dq (ν q + ω q )/2 − E q=q 0 .
Therefore, G(q) near q = q 0 can be written as which means that integral (122) becomes a Gaussian integral (with a complex coefficient A).In this way we arrive at In Fig. 8 we compare the exact result (126) with a numerical calculation, where we have used m = 2 , γ = 1 with a finite cutoff a = 1/20.The blue line presents the numerical result of the real part of (122) and the red line shows the real part of the analytic result (126), which is applicable above the energy E = E Uth = 1.632.We thus also see here an excellent agreement of the numerical result with the exact result.
FIG. 7. The shaded area is the positive Re G(q) region on the complex plane of q.The black point is q 0 (zero of (ν q + ω q )/2 − E).The orange line (contour) starts from the origin, goes through q 0 and approaches +∞, without leaving the shaded area.

V. CONCLUSION
In this paper, we have examined whether complex ghosts are truly not created by collisions of (positive norm) physical particles, as claimed by Lee and Wick [8][9][10].This problem is of very general interest because, if their claim is true, all theories can in principle be made renormalizable or even finite without violating unitarity, simply by adding higher derivative regulators.More importantly, quadratic gravity theory which comes with a natural built-in regulator, becomes a viable perturbatively renormalizable theory of gravity.To the author's knowledge, no clear disproof nor sound proof has been given for physical unitarity itself, aside from some severe criticisms on the Lorentz invariance [12,13] and causality [8,11,17,18].
Our answer to this question is: Complex ghosts can in fact be created with finite (nonzero) probability by the collisions of physical particles alone, thus implying that physical unitarity is violated.This is a clear and inevitable conclusion as far as one works faithfully within the QFT framework.It should also be emphasized that there is no room for inequivalent changes of the integration contours in the complex plane of the energy-momentum variables in Feynman diagrams.In other words, if one finds a clever method to change the integrations in such a way as to satisfy unitarity, it is no longer a quantum field theory.Furthermore, if one leaves the operator formalism of QFT, it would be very difficult to prove the unitarity in closed form.We were able to avoid the cumbersome (and sometimes difficult) considerations of the integration contours of the energy variables by using the 3d momentum form for the propagators, which amounts to a great simplification attained at the sacrifice of apparent covariance.We note, however, that it may be desirable to develop a more systematic computation technique by using the 4d expression propagators.
Although the energy conservation becomes slightly obscured by the appearance of a complex mass squared, we have found a clear threshold energy for the creation of the complex ghost, that is, the lower threshold given by E Lth = Re √ M 2 −Im √ M 2 .At energies below this threshold, there occurs no ghost production, meaning that unitarity and renormalizability hold.Consequently, we have a consistent effective QFT for E < E Lth .
The complex delta function has played a central role in this complex ghost theory, and its property as a distribution is very natural and mathematically consistent.Moreover, it can be used whenever states with complex energy eigenvalue appear, irrespective of the norm.Indeed, Nakanishi [26,27] introduced this complex delta function to construct a complex energy eigenstate corresponding to the usual (i.e., positive norm) unstable particle.It may thus be interesting to develop a QFT of unstable particles and compare it with the present complex ghost theory in order to answer the question: What is the essential difference between an unstable particle and a complex ghost?neither the anti-ghost pole ∓(ω * q − iε) nor the ghost pole ±(ω q + iε) by an infinitesimal amount iε.Therefore, we have to take care of the pole contribution only for the 'photon' case φ a = A in the second line integral in Eq. (A2): We can deform the contour R for the 'photon' propagator to for the terms containing δ c (P 0 − q 0 + ω q ′ ) and δ c (P 0 − q 0 − ω q ′ ), respectively.

(A7)
Here we have used the fact that the δ c (P 0 − q 0 ± ω q ′ ) reduces to the usual Dirac delta δ(P 0 − Re q 0 ± Re ω q ′ ) for q 0 on the shifted contour R(±ω q ), respectively.In going to the third equality, we have decomposed the preceding quantities into four partial fractions by (X 2 − Y 2 ) −1 = (2X) −1 ± (X ± Y ) −1 and recombined them into two terms.In the last line, we have kept iε only in the denominator which may become 0.

5 .
FIG.5.The unshaded region is the negative G R (q) region, i.e., where G R (q) < 0, and thus gives the support of the complex delta function δ c (E − ω q ).The upper and lower boundary curves of the negative G R (q) region are given by E = Re ω q ± Im ω q ≡ E ± (q).The curve in the middle of the region is E = Re ω q drawn for reference.This figure is drawn with parameter γ 2 /m 2 = 0.5.The lowest horizontal line (orange) is at lower threshold energy E = E − (q=0) ≡ E Lth = 0.786 and the upper horizontal line (green) is at the upper threshold energy E = E + (q=0) ≡ E Uth = 1.272.G R (q) = 0 has no (region I), one (region II) and two solutions (region III) of real q.

FIG. 6 .
FIG. 6.Comparison of the numerical calculation (blue) of Re I(E) with that of the exact result (red) (120), where m = 1 and γ = 1/ √ 5 with a finite cutoff a = 1/14 are used.The vertical lines separate the three energy regions, where the threshold values are: E Lth = 0.7862 , E Uth = 1.272.

FIG. 8 .
FIG. 8. Comparison of the numerical calculation (blue) of the real part of I ph (E) with that of the exact result (red) (126), where m = 2 and γ = 1 with a finite cutoff a = 1/20 are used.The vertical lines separate the three energy regions, where the threshold values are: E Lth = 1.384 , E Uth = 1.632.