Fakeons and Lee-Wick models

The “fakeon” is a fake degree of freedom, i.e. a degree of freedom that does not belong to the physical spectrum, but propagates inside the Feynman diagrams. Fakeons can be used to make higher-derivative theories unitary. Moreover, they help us clarify how the Lee-Wick models work. In this paper we study the fakeon models, that is to say the theories that contain fake and physical degrees of freedom. We formulate them by (nonanalytically) Wick rotating their Euclidean versions. We investigate the properties of arbitrary Feynman diagrams and, among other things, prove that the fakeon models are perturbatively unitary to all orders. If standard power counting constraints are fulfilled, the models are also renormalizable. The S matrix is regionwise analytic. The amplitudes can be continued from the Euclidean region to the other regions by means of an unambiguous, but nonanalytic, operation, called average continuation. We compute the average continuation of typical amplitudes in four, three and two dimensions and show that its predictions agree with those of the nonanalytic Wick rotation. By reconciling renormalizability and unitarity in higher-derivative theories, the fakeon models are good candidates to explain quantum gravity.


Introduction
The Lee-Wick (LW) models are a subclass of higher-derivative theories, where the free propagators contain complex conjugate pairs of extra poles, besides the poles corresponding to the physical degrees of freedom and the degrees of freedom due to the gauge fixing. The LW models are claimed to lead to a perturbatively unitarity S matrix [1][2][3] due to a certain compensation mechanism.
Various issues concerning the formulation of the LW theories remained open for a long time. For example, if they are defined as initially suggested by Lee [4], the models violate Lorentz invariance [5]. This problem is due to the incompleteness of the initial Lee-Wick prescription. Lee and Wick specified how to integrate on the loop energies, but did not provide a compatible prescription for the integrals on the loop space momenta.

JHEP02(2018)141
To overcome these difficulties, further prescriptions were supplemented later. For example, in ref. [3] a procedure of limit, which is known as CLOP prescription, was proposed to treat the critical situations where the LW poles pinch the integration paths on the complex energy planes. Lorentz invariance is recovered [6], but in some one-loop diagrams the CLOP prescription is ambiguous [7] and other ambiguities are expected at higher orders [3]. Moreover, it is unclear how to incorporate the CLOP prescription at the Lagrangian level or in the Feynman rules.
The problems were recently solved by reformulating the LW models by (nonanalytically) Wick rotating their Euclidean versions [7]. This procedure not only provides the correct prescription to integrate on the loop energies, which agrees with the Lee-Wick one, but also provides the natural companion prescription to integrate on the loop space momenta.
Briefly, the Lee-Wick integral on the loop energies includes complex values, so an integral on real values of the loop space momenta is not compatible with Lorentz invariance. However, if the integration domain on the loop space momenta is deformed in a suitable way to include complex values, Lorentz invariance is recovered.
It turns out that the Wick rotation is analytic only in a region of the space P of the (complexified) external momenta, the region that contains the purely imaginary energies. We call it main region and denote it by A 0 . The Wick rotation is nonanalytic elsewhere, due to the LW pinching [7]. In the end, the space P is divided into disjoint regions A i of analyticity. A loop integral gives an analytic function in each A i . The relations among the functions associated with different regions are unambiguous, but not analytic.
The domain deformation mentioned above is simple to formulate, but hard to implement practically. Fortunately, there exists a shortcut to get directly to the final result, which is simple and powerful. As said, the Wick rotation is analytic in A 0 . The obstacles that prevent the analytic continuation beyond A 0 are the LW thresholds, associated with LW poles that pinch the integration paths on the energies. The thresholds have the form p 2 =M 2 , where p is a linear combination of incoming momenta andM is a linear combination of (possibly complex) masses. A LW threshold can be analytically overcome in two independent ways. Neither of the two is separately compatible with unitarity and there is no way to choose between them. We show that the nonanalytic Wick rotation picks the arithmetic average of the two continuations, which we call average continuation. The final amplitudes are unitary, Lorentz invariant and analytic in every A i , i = 0, although not analytically related to the amplitudes evaluated in A 0 .
In this paper we study these issues in detail in arbitrary diagrams and show that the formulation of the LW models is consistent to all orders. We compute the average continuation of typical physical amplitudes in four, three and two spacetime dimensions and provide numerical checks that the average continuation and the nonanalytic Wick rotation give the same results. Moreover, we prove that the LW models are perturbatively unitary to all orders and show that their renormalization coincides with the one of their Euclidean versions. This property ensures that the locality of counterterms and the usual rules of power counting hold in every region A i .
The average continuation is an extremely powerful tool. It simplifies the computation of the amplitudes in the regions A i , i = 0. It eliminates the need of starting from the

JHEP02(2018)141
Euclidean space and performing the Wick rotation. It allows us to prove the perturbative unitarity to all orders in a relatively straightforward way. It gives an efficacious control on the renormalization.
In ref. [8] the perturbative unitarity of the LW models was proved at one loop. The generalization of the proof to all orders can be worked out by first deriving the so-called cutting equations [9][10][11][12] (which imply the unitarity equation SS † = 1), in the main region A 0 and then proving that they can be average-continued to the other regions A i . The final cutting equations have the expected, unitary form and propagate only the physical degrees of freedom. We actually need to work with generalized versions of the equations, which are proved starting from the algebraic cutting equations (ACE) of ref. [12], a set of polynomial identities associated with Feynman diagrams which are particularly fit to perform the average continuation from A 0 to A i .
We recall that the cutting equations imply SS † = 1 straightforwardly in the models involving just scalar fields and fermions. In gauge theories [13,14] and gravity [11], they imply a pseudounitarity equation, which turns into the unitarity equation after proving that the temporal and longitudinal components of the gauge fields are compensated by the Faddeev-Popov ghosts.
It is important to stress that not all the higher-derivative theories fall in the Lee-Wick class. For example, the Lee-Wick models of quantum gravity are typically superrenormalizable. The reason is that the LW poles must come in complex conjugate pairs, which requires many higher derivatives. With fewer higher derivatives we may build a strictly renormalizable theory [15][16][17][18], but then the free propagators have ghost poles with real squared masses. In ref. [19] it was shown that it is possible to double such poles by means of a new quantization prescription and treat them as LW poles associated with a fictitious LW scale E that is sent to zero at the very end. This leads to the introduction of the notion of fake degree of freedom, or "fakeon". Once a pole is doubled according to this prescription, it can be consistently dropped from the physical spectrum. Turning ghosts into fakeons allows us to make the higher-derivative theories unitary.
The notion of fakeon generalizes the ideas of Lee and Wick and actually clarifies their crucial properties. For example, the nonanalyticity of the S matrix due to the LW pinching can be seen as associated with a fakeon of a finite LW scale E = M . For this reason, the LW models are particular "fakeon models", by which we mean models with physical degrees of freedom and fakeons. The results of this paper, such as the proof of perturbative unitarity to all orders, hold in all the fakeon models.
We recall that the LW models are also investigated for their possible phenomenological implications, for example in QED [2], the standard model [20][21][22][23] and grand unified theories [24,25], besides quantum gravity [19,[26][27][28][29]. The results of this paper and refs. [7,8,19] raise the fakeon models to the status of consistent fundamental theories, since the theoretical problems that could justify a certain skepticism around them are now overcome. In particular, we have viable candidates to explain quantum gravity within quantum field theory. Among the various possibilities, a unique one is strictly renormalizable [19].
The paper is organized as follows. In sections 2 and 3 we recall the formulation of the Lee-Wick models as nonanalytically Wick rotated Euclidean theories and investigate JHEP02(2018)141 their main properties in arbitrary Feynman diagrams. In particular, in section 2 we study the LW pinching, while in section 3 we study the domain deformation. In section 4 we define the average continuation of an analytic function and analyse its properties. We also define the difference continuation, which is useful for the cutting equations. In section 5 we study the average continuation of typical amplitudes in various dimensions and numerically compare the results with those of the nonanalytic Wick rotation. In section 6 we recall the definition of fakeon and it main properties. In section 7 we prove the perturbative unitarity of the fakeon models to all orders. In section 8 we show that the counterterms of the fakeon models are the same as those of their Euclidean versions. Section 9 contains the conclusions.

Lee-Wick models
In this section we study the Lee-Wick models by nonanalytically Wick rotating their Euclidean versions. The arguments hold to all orders in spacetime dimensions D greater than or equal to two, in local quantum field theories whose free propagators have poles that are located symmetrically with respect to the real axis of the complex energy plane, with squared masses that have nonnegative real parts. The poles located on the real axis are called standard poles and the other ones are called LW poles. The standard poles are physical if they have positive residues.
Observe that derivative vertices and propagators with nontrivial numerators do not change the analysis that follows. What matters in a loop integral are the singularities of its integrand, i.e. the denominators of the propagators.
Before plunging into the nonanalytic Wick rotation, let us stress why alternative approaches to the formulation of higher-derivative theories are not viable. Letting aside ad hoc prescriptions such as the CLOP one, which cannot be incorporated at the level of the Feynman rules and lead to ambiguous results, a natural formulation that may come to mind is the Minkowski one, where the loop energies are integrated on their natural, real values. Recently, it has been shown that the Minkowski formulation generates nonlocal, non-Hermitian divergences that cannot be removed by any standard procedures [30]. In the few cases where the locality of counterterms is not violated, the amplitudes are not consistent with perturbative unitarity [8]. These observations lead us to conclude that the Minkowski formulation is not the right one. The only chance to define the higher-derivative models consistently is the Wick rotation of their Euclidean versions.
The simplest example of LW propagator is where M and µ are real mass scales. The poles of this propagator are shown in figure 1. The standard poles are encircled and read p 0 = ±ω (p), where ω (p) = p 2 + m 2 − i and p = (p 0 , p). The LW poles are not encircled and read p 0 = ±Ω + (p) and where Ω ± (p) = p 2 + M 2 ± and M ± = µ 2 ± iM 2 . We call the pairs of poles Ω ± and −Ω ± Lee-Wick pairs. Note that the Minkowski and Euclidean versions of the theories are JHEP02(2018)141  not equivalent, since the free propagators have poles in the first and third quadrants of the complex plane.
Following ref. [7], the loop integrals are defined starting from the Euclidean version of the theory. In the case of the tadpole diagram, the Wick rotation leads to the integration path shown in figure 1. We see that the poles that are located to the right (resp. left) of the imaginary axis are below (above) the integration path.
The bubble diagram which involves the product of two propagators, better illustrates the general case. There, the Wick rotation leads to integration paths of the form shown in the left picture of figure 2. The thick crosses denote the poles of the propagator S(k − p, m 2 ), which depend on p. The other crosses denote the poles of S(k, m 1 ), which are fixed. The general rule, which holds for arbitrary diagrams, is that the right (resp. left) poles of a propagator -i.e. those whose energies have positive (negative) real parts at zero external momenta -are located below (above) the integration path.
When we wary p, a LW pole of S(k − p, m 2 ) can approach a LW pole of S(k, m 1 ) from the opposite side of the integration path. When the two come to coincide, we have a Lee-Wick pinching. The standard poles can give the usual pinching, which we call standard pinching. Similarly, a mixed LW pinching involves a LW pole and a standard pole.
The condition for having a LW pinching is a system of two pole conditions. For example, the right picture of figure 2 describes the simultaneous pinching of the poles of two LW pairs. The conditions for the top pinching are while the conditions for the bottom pinching are their complex conjugates (with the understanding that the conjugation does not act on the momenta). Solving (2.3) for k 0 , we obtain Varying k in R 3 with p real and fixed, the solutions of this equation fill the region enclosed inside the curve γ of figure 3. Other LW pinchings occur for and fill the regions enclosed inside the other two curves of figure 3. Finally, we have the regions obtained by reflecting (2.4) and (2.5) with respect to the imaginary axis. Summarizing, the complex plane is divided into certain regions, which we denote bỹ A i . The curve γ is the boundary of the regionÃ P that intersects the positive real axis. The region that contains the imaginary axis is denoted byÃ 0 .
The regionsÃ i are not Lorentz invariant, which is the reason why they are not the final analytic regions A i . For example, the threshold of the LW pinching given by eq. (2.4) is the point P with p 2 = 2µ 2 + 2 µ 4 + M 4 ≡ 2M 2 LW , as we prove below. However, the intersection between the curve γ and the real axis is not P , but a different point P . It is useful to introduce two functions Then the point P has energy This relation cannot be expressed as a Lorentz invariant threshold condition of the form For a while, we focus on real external momenta p, which are the ones of physical interest. Note that (2.7) satisfies 4µ 2 p 2 2M 2 LW , where the equalities holds for p 2 = ∞ and p = 0, respectively.
We define the Euclidean region as the strip |Re[p 0 ]| < |p|, which contains the imaginary axis. It is easy to check that the LW pinching conditions do not admit solutions there. Indeed, formulas (2.4) and (2.5) show that when a LW pinching occurs, the minimum of |Re[p 0 ]| is the right-hand side of (2.7), which is greater than or equal to p 2 + 4µ 2 . In particular, the Euclidean region is a subregion ofÃ 0 .
We define the loop integral B(p) as follows. First, we integrate on the loop energy k 0 by means of the residue theorem. Then, we concentrate on the Euclidean region and integrate the loop space momentum k on its natural domain R 3 . Since no LW pinching occurs, the result is analytic (and Lorentz invariant) but for the branch cuts associated with the standard pinching.
Next, we ask ourselves if we can analytically extend the result away from the Euclidean region. Focusing on the real axis, we find no obstacle for p 2 < 4µ 2 , because all such points are below P . We can also reach values p 2 4µ 2 , as long as we restrict the Lorentz frame to the subset where the LW pinching does not occur for any k ∈ R 3 . The good frames are those that have energies p 0 smaller than the energy of P . By formula (2.7), this condition can be written as which admits solutions if and only if p 2 < 2M 2 LW (with p 2 4µ 2 ). In the end, for p 2 < 2M 2 LW , there is always an open subset L of Lorentz frames where no LW pinching occurs and we can evaluate the loop integral by integrating k on R 3 . The result is the analytic continuation of the function obtained in the Euclidean region. Since it does not depend on the Lorentz frame, it can be straightforwardly extended from L to the whole space of Lorentz frames.
We have thus proved that the true LW threshold is the point P of figure 3, beyond which the LW pinching is inevitable and the regionÃ 0 cannot be extended further. The region A 0 , which is the maximal extension ofÃ 0 , stops at P .
The true challenge of the Lee-Wick models is to overcome the LW threshold P . To make a step forward towards the solution of this problem, we generalize the calculation just described as follows. So far, we have calculated the loop integral in a specific subset L of Lorentz frames, for 4µ 2 < p 2 < 2M 2 LW , because we wanted to be able to integrate k on R 3 . Then, we extended the result to all the Lorentz frames by Lorentz invariance. If we want to make the calculation for 4µ 2 < p 2 < 2M 2 LW directly in an arbitrary Lorentz frame,

JHEP02(2018)141
we must deform the k integration domain D k to ensure that the LW pinching does not occur for any p 2 . For example, if O P denotes the portion of the real axis with p 2 2M 2 LW , p 0 > 0, we can choose a deformation that squeezes the regionÃ P onto O P (see the next section for details). Observe that O P is Lorentz invariant.
The good news is that the domain deformation just mentioned allows us to work out the loop integral even beyond the LW threshold P . In that case, we have to proceed as follows. LetÃ def P denote the deformed regionÃ P , before it is squeezed onto O P . Let D def k denote the k integration domain associated withÃ def P . We go insideÃ def P and evaluate the loop integral B(p) there. Since the condition (2.4) is complex, it can be split into two real conditions x = y = 0 for suitable functions x and y of k. Changing variables, in D 3 the singularity has the form dxdy x + iy , (2.9) which is integrable. In D = 2 there is no singularity, because the pinching just occurs at the boundaries γ, γ def of the regionsÃ P ,Ã def P . We view the result of the calculation iñ A def P as a function of the k integration domain D def k . When we finalize the deformation that takesÃ def P to O P , we obtain the value of the loop integral on the real axis above the LW threshold P .
At the end, we can take A P = O P . Alternatively, we can analytically continue the result found in O P to a neighborhood of O P and take that neighborhood as the final region A P (reducing A 0 correspondingly).
Before the squeezing ofÃ P to O P , the result of the loop integral inÃ P is neither analytic nor Lorentz invariant, in dimensions greater than or equal to three. On the other hand, two dimensions are exceptional, because in D = 2 the LW pinching occurs only at the boundaries of the regionsÃ i , i = 0, but not inside. Consequently, the loop integral is both Lorentz invariant and analytic in (a neighborhood of) O P , even before making the domain deformation. We check these properties explicitly in the examples of section 5. In the next section we explain in detail how the domain deformation works in arbitrary D 2.
So far, we have focused on the LW thresholds that are located on the real axis. Similar arguments hold for the other LW pinchings (2.5), whose thresholds are the points of minimum Re[p 0 ] of the corresponding regionsÃ i . It is easy to check that such points have Re[p 0 ] = 2η + (p 2 /4 + µ 2 ) and Im[p 0 ] = ±2η − (p 2 /4 + µ 2 ), so the thresholds are p 2 = 4µ 2 ± 4iM 2 . When p → 0 the corresponding regionsÃ i squeeze onto curves with endpoints at the thresholds. The calculations beyond such thresholds are performed with a procedure analogous to the one described above: first, we evaluate the loop integral inside a regionÃ i ; then, we deform the k integration domain tillÃ i gets squeezed onto a curve; finally, we take such a curve as the final region A i , or enlarge it to some neighborhood of it by analytic continuing the result found in it.
More complicated one-loop diagrams can be studied similarly. As an example, consider the box diagram shown in the left picture of figure 4. We assume that the propagators have the same masses m, µ and M , for simplicity. The pinchings may occur when two, three or four propagators have simultaneous pole singularities. Decomposing into partial fractions, the integrand can be written as a sum of terms , (2.10) where z denotes the loop energy k 0 , σ i = ±1 and each a, b, c, d is a frequency ω or Ω ± plus a linear combination of incoming external energies. The poles with σ i = 1 lie one side of the z integration path, while the poles with σ i = −1 lie on the other side. If all the σ i are equal, the residue theorem gives zero. If σ 1 = 1, σ 2 = σ 3 = σ 4 = −1, the residue theorem give a result proportional to 1 (a + b)(a + c)(a + d) . (2.11) If σ 1 = σ 3 = 1, σ 2 = σ 4 = −1, the residue theorem gives . (2.12) Each singularity of (2.11) and (2.12) has the form (2.4) or (2.5). The other cases are permutations of the ones just described. Note that the frequencies are always summed with positive relative signs.
In the end, we only have situations that are analogous to those already met in the case of the bubble diagram. The LW thresholds are where p i , i = 1, 2, 3 denote the incoming momenta shown in the picture. The evaluation of the loop integral proceeds as before. We first compute it in the Euclidean region, where no LW pinching occurs, by integrating the loop space momentum k on R 3 . Then we extend the result by analytic continuation toÃ 0 . Third, we maximize the regionÃ 0 , again by analytic continuation, which identifies the region A 0 . Beyond A 0 we find obstacles, given by the LW thresholds. We overcome those obstacles by going inside the regionsÃ i , i = 0, and then deforming the k integration domain to squeeze those regions into curves. At the end, we may define the regions A i as neighborhoods of those curves. We arrange each A i so as to make it Lorentz invariant for real external momenta.

LW pinching beyond one loop
Before considering an arbitrary multiloop diagram, we begin with the chestnut diagram shown in the right picture of figure 4. The propagators 1 and 2 depend on one loop momentum, which we call k. The integration path over k 0 gets pinched when two poles come to coincide from opposite sides. This gives relations of the form whereω i can stand for ω or Ω ± and p is an external momentum. Integrating over k 0 by means of the residue theorem, we remain with a single pole, which occurs for This condition is analogous to (2.4) and (2.5). Now, let us consider the propagators 1, 3 and 4. They depend on two loop momenta, k 1 and k 2 , which we assign to the legs 1 and 4. Their simultaneous singularities give pinching conditions of the form where p is a sum of incoming external momenta. The signs in front of the frequencies ensure that the first and third pole lie on opposite sides with respect to the k 0 1 integration path, while the second and third pole lie on opposite sides with respect to the k 0 2 integration path. The integrals over k 0 1 and k 0 2 eliminate the first two conditions and turn the third one into p 0 =ω 1 (k 1 ) +ω 2 (k 2 ) +ω 3 (k 1 + k 2 − p). (2.14) Now, let us consider the contribution where z 1 = k 0 1 , z 2 = k 0 2 and a, b, c and d are defined as before and associated with the legs 1, 2, 3, and 4, respectively. A pinching can occur, since a lies on one side of the z 1 integration path, with b, c on the other side of it, and at the same time c and d lie on opposite sides of the z 2 integration path. The residue theorem gives a result proportional to .
The denominator vanishes in three situations, two minimal and one nonminimal. The minimal condition a + b = 0 has the form (2.13). The minimal condition a + c + d = 0 has the form (2.14). The nonminimal condition is the system made of the two. The calculation can proceed as in the one loop case, the only difference being that at some point we have to deform the integration domains of both loop space momenta. The other contributions to the chestnut diagram can be treated similarly.
The arguments just given can be generalized to diagrams with arbitrary numbers of loops. The minimal configuration of pole singularities which may give a pinching occurs JHEP02(2018)141 when the number n of propagators that have simultaneous pole singularities is equal to the number of loop momenta they depend on, plus one. If we parametrize the loop momenta in a convenient way, the first n − 1 conditions read k 0 i =ω i (k i ), i = 1, . . . n − 1. After integrating on the loop energies k 0 i by means of the residue theorem, the last condition becomes and p is again a sum of incoming external momenta. This is the minimal pinching condition, with a convenient parametrization for the momenta. More generally, the k i may be independent linear combinations of the loop momenta (with coefficients ±1) plus linear combinations of the external momenta.
The most general configuration of pole singularities arises as a superposition of minimal configurations (plus configurations of singularities that give no pinching, which we can ignore). Then, the most general pinching condition is just a system made of minimal conditions. For this reason, it is sufficient to study the minimal condition, in the parametrization (2.16).
We may have a pure LW pinching, where only LW poles are involved, a mixed LW pinching, where both LW and standard poles are involved, and a standard pinching, where only standard poles are involved.
An important fact is that the signs in front of the frequencies that appear on the right-hand side of (2.16) are always positive. The reason is that the pinching just occurs between right and left poles of different propagators, the right ones being placed below the integration path on the loop energy and the left ones being placed above it. There is no pinching between two right poles or two left poles (which would generate minus signs in front of the frequencies), because they are located on the same side of the integration path.
The threshold associated with the pinching condition (2.15)-(2.16) is (2.17) This formula is a straightforward generalization of the one that holds in the standard case, but must be proved anew, because the LW pinching involves unusual features, such as the extended regionsÃ i that violate Lorentz invariance in some intermediate steps.
Specifically, the thresholds are found by means of a two-step procedure: first we minimize Re[p 0 ] in k i and then we maximize Re[p 2 ] in p. Referring to the analysis made at one loop for thresholds on the real axis, the first step corresponds to identifying the point P of figure 3 and the second step corresponds to deforming P into P . Now we prove that this procedure does give formula (2.17).

JHEP02(2018)141
Let us first consider the case where only LW poles are involved, i.e. n + frequenciesω i are equal to Ω + and n − frequenciesω i are equal to Ω − , with n = n + + n − . We have where η + is defined in formula (2.6) and Minimizing Re[p 0 ] in k i , we obtain k i = p/n for every i, which gives The maximum of Re[p 2 ] in p is its value for p = 0, which gives the thresholds The result agrees with (2.17), since Ω ± (0) = η + (µ 2 ) ± iη − (µ 2 ). The thresholds on the real axis are those with n + = n − . Observe that no LW pinching occurs in the Euclidean region |Re[p 0 ]| < |p|. Indeed, using formula (2.19a) we find that wherever a LW pinching occurs the inequalities Next, let us consider the mixed LW pinching, where both standard poles and LW poles are present. We assume that µ and M are the same everywhere, but the standard masses are generic. We separate the last standard pole, with mass m, from the other ones, with masses m j and loop space momenta q j . Then, we get the condition (2.15) with (2.20) Here we have defined ω(p, m) = p 2 + m 2 and First, we minimize Re[p 0 ] in q, which is straightforward. Indeed, translating Re[p 0 ] by a constant, this operation just gives the threshold of the standard pinching. We thus find

JHEP02(2018)141
and Now we minimize Re[p 0 ] in k i , which gives k i = pα(p) ≡ s for every i, for some function α of p. It is convenient to express everything in terms of s rather than p. We find where η + (x) = dη + (x)/dx. Unless specified differently, here and below the arguments of η + , η − and their derivatives are s 2 + µ 2 . It is easy to check that the argument of the square root in (2.22) is always positive. Formula (2.21) gives q j = m j (p − ns)/m tot . Using (2.22) inside D pinch = 0, we get At this point, we maximize Re[p 2 ] in p. We can actually maximize it in s, since dp 2 /ds 2 is always positive. It is easy to show that the right-hand side of (2.23b) is a monotonically decreasing function of s 2 , so the maximum of Re[p 2 ] coincides with its value at s = 0, which gives the threshold in agreement with (2.17). Again, no LW pinching occurs in the Euclidean region |Re[p 0 ]| < |p|. Indeed, for arbitrary k i and q j , the LW pinching conditions D pinch = 0 imply is the momentum p that minimizes Re[p 0 ] in k i and q j , encoded in formulas (2.23a) and (2.23b).
Consider a Feynman diagram G with n + 1 external legs. Let p 1 , · · · , p n denote the incoming momenta of n external legs. The thresholds read where I is a subset of indices of the incoming momenta andM is positive sum of ordinary masses m and LW masses M ± . Note that the incoming momentum of the (n+1)th external leg is so whenever a sum of incoming momenta includes p n+1 it can be written as minus a sum of p i . Since the overall sign is immaterial for the left-hand side of (2.25), we can always write the thresholds as in that formula. The number of thresholds (2.25) and regionsÃ i of each loop integral is finite. If the masses m j are nonvanishing and finitely many, the number of thresholds of an amplitude is finite within any compact energy range, even after summing the loop corrections to all orders. That number becomes infinite when some masses m j vanish. This is the known problem of the infrared divergences, which is dealt with by means of resummation techniques [31][32][33][34].
Strictly speaking, the sum m tot of standard masses in formula (2.24) should be equipped with a small negative imaginary part, coming from the width of the propagator (2.1). In several calculations, as well as the proof of perturbative unitarity of section 7, it is necessary to work at = 0. Then the thresholds (2.24) with n + = n − are not exactly on the real axis for m tot = 0, but a bit displaced from it. As before, when LW poles are involved, the conditions (2.15) identify extended regionsÃ i , i = 0. Since is supposed to be small, while M is finite, the regionsÃ i always intersect the real axis in a segment, when n + = n − . A typical situation is shown in figure 5, where P 1 and P 3 are standard thresholds, while P 2 and P 4 are LW thresholds. For convenience, we have drawn the branch cuts ending at the standard thresholds so that they do not intersect the regionsÃ i , i = 0.
A loop integral I is first evaluated in the Euclidean region, by integrating on the natural real domain R 3(n+r−1) of the loop space momenta k i and q j . Then the result is extended by analytic continuation toÃ 0 and A 0 . Above the LW thresholds, the integration domain D k,q on k i and q j is deformed from R 3(n+r−1) till the regionsÃ i , i = 0, squeeze onto Lorentz invariant surfaces L i . The calculation of I is performed inside each deformed A i , i = 0, before finalizing the squeezing. Once the squeezing is finalized, the results found in the surfaces L i are extended to neighborhoods of them by analytic continuation. Those neighborhoods can be taken as the regions A i , i = 0. For every threshold with n + = n − , the corresponding region A i is enlarged enough till it intersects the real axis in a segment, as in figure 5. Note that the singularities 1/D pinch associated with the LW pinchings have the form (2.9) and so are integrable.

The domain deformation
In the most general case, the deformation of the integration domain on the loop space momenta, required by the nonanalytic Wick rotation, is a rather involved process. However, its main features are relatively simple. In this section we illustrate them in detail, starting JHEP02(2018)141 from the bubble diagram in D = 2, then generalizing the arguments to arbitrary D and arbitrary diagrams.

Domain deformation in the bubble diagram
Consider the LW pinching condition (2.4) Let us keep p x fixed (and real) and view k ± x as functions of p 0 . If we move p 0 on figure 3 along lines parallel to the real axis, with Im[p 0 ] 0 and Im[p 0 ] 0, we obtain the pictures of figure 6 (where M = µ = p x = 1). In each picture, the trajectories are the functions k ± x (p) and the arrows point towards growing values of Re[p 0 ]. As long as Im[p 0 ] = 0, the trajectories do not intersect each other. If we take the limit Im[p 0 ] → 0, we obtain figure 7, where the points a i with the same index i correspond to solutions k x with the same value of p 0 .
The natural k x integration domain is the k x real axis. In this discussion we denote it by D 1 . Let us follow the solutions of figure 7 and see how the integration domain must be deformed to have analyticity. Referring to figure 3, we start from the segment of the p 0 real axis that is located below P . A typical point there is sent into the two points a 1 JHEP02(2018)141 Figure 8. Basic domain deformation. of figure 7, which are located on opposite sides of the domain D 1 . When p 0 increases, one trajectory k ± x intersects D 1 (which happends when p 0 reaches the point P ). The segment of the p 0 real axis contained between P and P is represented by the points a 2 in the k x complex plane, which are located on the same side of D 1 . The loop integral, as a function of p 0 , can be analytically extended beyond P by deforming D 1 into some new domain D 2 that looks like the one shown in the first picture of figure 8, so that the points a 2 are left on opposite sides.
When p 0 continues to increase and reaches the point P , the two trajectories hit each other. There, it is impossible to keep the solutions k ± x on opposite sides of the k x integration domain. This means that the loop integral cannot be analytically extended beyond P by moving p 0 along the real axis. The point P is the sole and true case where the pinching cannot be avoided. It is obtained by setting the argument of the square root of (3.1) to zero, which gives the LW threshold p 2 = 2M 2 LW , in agreement with the results of the previous section.
Larger real values of p 0 take us into the portion O P of the real axis above P , which is represented by the points a 3 of figure 7. There are two types D 3 and D 3 of deformed domains that leave those points on opposite sides, as shown in the second and third pictures of figure 8. The two possibilities correspond to reaching O P by giving p 0 a small positive, or small negative, imaginary part. Indeed, we know from figure 6 that the analytic continuation finds no obstacles in those cases, because the k x trajectories never intersect each other.
In the end, we have two analytic continuations fromÃ 0 to O P , one obtained by circumventing P from the half plane Im[p 0 ] > 0 and the other one obtained by circumventing P from the half plane Im[p 0 ] < 0. We will see in section 4 that the result of the loop integral above P is the arithmetic average of the two (average continuation).
Finally, the regionÃ P can be completely squeezed onto O P by deforming D 1 into the domain D P made of the curve that crosses the points a 3 of figure 7. Indeed, figure 9 shows that D P always leaves the solutions k ± x on the same side, no matter how small |Im[p 0 ]| is taken.
The arguments can be easily extended to arbitrary dimensions D greater than two. Assume that the external space momentum p is directed along the x direction. Writing p = (p 0 , p x , 0) and k = (k 0 , k x , k ⊥ ), it is easy to check that the conditions (2.4) and (2.5) in D > 2 are obtained from those in D = 2 by means of the replacement µ 2 → k 2 ⊥ + µ 2 ≡μ 2 . Then it is apparent that to squeeze the regionÃ P onto O P we do not need to deform k ⊥ to complex values, since it is enough to deform the k x integration domain as explained above, for everyμ 2 .

JHEP02(2018)141
To summarize, the equations (2.4) and (2.5) tell us when the integration path on the loop energy gets pinched. However, in most cases the pinching is eventually avoided by deforming the integration domain on the loop space momenta. The pinching is inevitable only at the LW thresholds. Since the LW thresholds are Lorentz invariant, Lorentz invariance is never truly violated. Moreover, the regionsÃ i , i = 0, can be deformed and squeezed at will. The regions located above the LW thresholds can be reached analytically from the regions located below the LW thresholds in two independent ways.
It should also be noted that everything we have said so far equally applies to the standard thresholds and actually offers a new approach to investigate their properties. In the limit M → 0 the solutions (3.1) become In particular, we can appreciate why the thresholds are the only points of true pinching, while the points lying on the branch cuts are not. Indeed, the branch cuts can be displaced at will by deforming the integration domains on the loop space momenta.

Domain deformation in more complicated diagrams
Now we study the domain deformation in the diagrams with more loops and/or more independent external momenta.
If we have a single threshold, the analysis of the previous subsection can be repeated with straightforward modifications. A unique combination p of external momenta is involved. If the pinching conditions involve a unique loop momentum k, the analysis is exactly the same as before. If they involve more than one loop momenta, we simply have more freedom to perform the domain deformation.
Thus, we can concentrate on the case of multiple LW thresholds. We begin from two LW thresholds involving the same combination p of external momenta. We denote them by P 1 : p 2 =M 2 1 , P 2 : p 2 =M 2 2 , etc. Let us assume the worst scenario (which reduces the freedom to make the deformation to a minimum), where there is a unique loop momentum k. As before, we choose p = (p 0 , p x , 0), k = (k 0 , k x , k ⊥ ). Consider the condition (2.4)  It may be helpful to see what happens with the help of a sort of animation. Then we see that, say, the points a 1 , b 1 lying on the trajectories that approach the threshold P 1 arrive first, while the points a 2 , b 2 lying on the trajectories that approach P 2 arrive later, as shown in figure 11.

JHEP02(2018)141
In figures 10 and 11 the symbols D P i , i = 1, 2, denote the k x integration domains that would squeezeÃ P i onto the real axis if the threshold were only P i . In the presence of both thresholds, we deform the k x integration domain into a "dynamic" domain D dyn (i.e. a function of p 0 ) as follows. At a first stage, when a 1 , b 1 approach P 1 and a 2 , b 2 are far away (in a neighborhood of the vertical Re[k x ] = p x /2), D dyn can be taken to be D P 1 . At a second stage, when a 1 , b 1 are far away in a neighborhood of D P 1 and a 2 , b 2 are approaching P 2 , we gradually deform D P 1 into D P 2 , starting from the vertical line towards the sides, as shown in figure 11. What said about the trajectories displayed in figures 10 and 11 can be repeated for the mirror trajectories obtained by reflection with respect to the vertical line, which correspond to the case Im[p 0 ] 0. Deforming the k x integration domain into D dyn as explained, no pinching ever occurs in the complex k x plane as long as |Im[p 0 ]| is sufficiently small and nonvanishing. This means that the domain deformation can be finalized as expected, till the regionÃ P squeezes completely onto the real axis of the complex p 0 plane.

JHEP02(2018)141
If the condition (2.4) with µ → µ 1 , M → M 1 is combined with the complex conjugate of the condition (2.4) (where the conjugation does not act on the momenta) with µ → µ 2 , M → M 2 , then P 2 and the trajectories approaching P 2 are reflected with respect to the real axis and with respect to the vertical line Re[k x ] = p x /2. The conclusions reached above can easily be extended to this case.
It can also be seen that the branch points due to the square roots involved in the expressions (2.4), (2.16) and (2.20) of D pinch are located away from the real axis of the k x complex plane (if µ 2 1 + k 2 ⊥ + M 2 1 and µ 2 2 + k 2 ⊥ + M 2 2 are nonvanishing, which we may assume here). Thus, if we choose |p x | large enough their branch cuts do not intersect the trajectories and domains described so far. Now we consider the case of two LW thresholds P 1 and P 2 that depend on different combinations p and q of external momenta, respectively. Again, we assume the worst scenario for the loop momenta, which is when only one of them is involved. This situation occurs, for example, in the triangle diagram. In D = 2 we have a picture such as the one of figure 12.
We see that the two domains D P 1 and D P 2 may intersect in a point I A , which is another true pinching. This kind of pinching also occurs in ordinary models, where it gives the socalled anomalous threshold [35]. In two dimensions the anomalous threshold of the triangle diagram is just a pole, but in higher dimensions it is a branch point. Other intersections that may give anomalous thresholds are those between D P 1 and the vertical line crossing P 2 , as well as the intersection between D P 2 and the vertical line crossing P 1 .
Anomalous thresholds are known to appear in the diagrams that involve more than one independent external momentum and have been studied at length in the triangle and box diagrams. Basically, any time there are two external momenta p and q, or more, singularities of the form ∼ 1/f (p 2 , q 2 , p · q) may appear, where f (p 2 , q 2 , p · q) is a nontrivial function of the invariants that can be built with them. Anomalous thresholds are associated JHEP02(2018)141 Figure 13. Domain deformation in the presence of multiple thresholds with different external momenta p and q. The grey dots are the points k ± x (p) and k ± x (q).
with cuts that split the diagram in more than two parts. It is known that they do not conflict with unitarity in ordinary models. We will see that this property extends to the Lee-Wick models. Ultimately, anomalous thresholds are sources of further complications, but do not pose new conceptual challenges. The dynamical squeezing can be achieved as follows. Consider the union D P 1 ∪ D P 2 and write it as D + ∪ D − , where D + (resp. D − ) is made of the superior (inferior) portions of D P 1 and D P 2 up to I A . Start from the domain D + . Consider the four trajectories k ± x (p) and k ± x (q) and take energies p 0 and q 0 that make them stay in neighborhoods of D + . Let p 0 and q 0 grow till the trajectories approach I A . If the trajectories k ± x (q) arrive first and the trajectories k ± x (p) arrive second, gradually deform D + into the domain show in figure 13. If k ± x (p) arrive first and k ± x (q) arrive second, take a domain deformation that is symmetric to the one of figure 13 with respect to the vertical line crossing I A . The two possibilities correspond to the two ways of circumventing the anomalous threshold I A . When p 0 and q 0 grow more, it is enough to stretch the deformations just described.
The arguments given so far easily extend to D > 2 and are exhaustive enough to understand what happens in the most general case.

Average continuation and difference continuation
When we start from the Euclidean version of the theory and perform the nonanalytic Wick rotation, we must deform the integration domain on the loop space momenta to overcome the LW thresholds. The domain deformation, described in the previous section, is not easy to implement in general. Fortunately, there is a shortcut to avoid it, which is the average continuation.
In this section we formulate the average continuation and show that it solves the nonanalytic Wick rotation and actually makes it unnecessary. Precisely, the average continuation allows us to calculate the loop integrals everywhere starting from the Euclidean region, or the regionÃ 0 , without even entering the other regionsÃ i , i = 0. We also study the difference continuation, which is an elaboration of a rather familiar concept, but helps clarify the properties of the average continuation by comparison.
The average continuation and the difference continuation are two noticeable nonanalytic procedures to define a function of a complex variable z beyond a branch point P . The JHEP02(2018)141 Re[z] Figure 14. Average continuation.
average continuation has to do with fakeons and ultimately solves the Lee-Wick models. The difference continuation is at the root of the cut diagrams. For simplicity, let us assume that P is located at the origin z = 0. Let f (z) denote the function we want to continue, defined by choosing the branch cut to be the negative real axis.
Referring to figure 14, define two other functions, f + (z) and f − (z), by choosing their branch cuts on the positive and negative imaginary axes, respectively, i.e. z = iρ and z = −iρ, with ρ 0. The average-continued function f AV (z) is defined as the average of f + (z) and f − (z): The imaginary axis divides the complex plane into two disjoint regions. This means that f AV (z) is actually a collection of two analytic functions: a superior function f > (z) = f AV (z)| Re[z]>0 and an inferior function The difference continuation is instead Clearly, f d (z) = 0 in the half plane Re[z] > 0. Among the properties of the average and difference continuations, we mention that: (i) the inferior function f < (z) is uniquely determined by the superior function f > (z), albeit in a nonanalytic way; (ii) the superior function f > (z) may or may not be determined by the inferior function f < (z); (iii) the superior function cannot be analytically continued beyond P ; (iv) it may or may not be possible to analytically continue the inferior function beyond P .
(v) if g(z) is analytic or has a pole in P and h(z) ≡ f (vi) if f (z) is real on the positive real axis, then f AV (z) and f d (z) are, respectively, real and purely imaginary on the real axis.
In the case (vi), the value of f AV (z) on the negative real axis is equal to the real part of either analytic continuation of f (z) to that half axis. Then we write on the whole real axis.
More generally, if f (z) has more distinct branch points on the real axis, the average and difference continuations are defined by applying the rules listed above to each branch point at a time. For example, let us study the average continuation with two branch points P 1 and P 2 . We have the situation depicted in figure 15, which leads to three disjoint regions: the half plane     > (z) for the second step and apply the average continuation again, to overcome P 2 and go from A 2 to A 3 . So doing, we obtain the new inferior function f When the branch points coincide, we must first deform them to make them distinct (by varying the masses), then apply the procedure just described and, finally, take the limit that makes them coincide. For example, consider a diagram G made of two diagrams G 1 and G 2 with one vertex in common. The average-continued function G AV (z) associated with G must clearly be the product G 1AV (z)G 2AV (z) of the average-continued functions G 1AV (z) and G 2AV (z) associated with G 1 and G 2 . However, if G 1 and G 2 have coinciding branch points, it may be tricky to satisfy this property. Consider figure 15 again: if P 1 = P 2 , we miss the paths shown in the second and third figure, so we may obtain a wrong continuation. For example, if f (z) = √ z is the superior function associated with G 1 = G 2 , then f 2 (z) = z is the superior function associated with G. However, z has no JHEP02(2018)141 branch point at all, rather than having two coinciding branch points, so it cannot give the correct result. Instead, if we replace f 2 (z) with √ z √ z − a, with a > 0, perform the average continuation and let a tend to zero at the end, we get the correct result. The outcome is independent of the deformation. Indeed, if we exchange the points P 1 and P 2 in figure 15, we simply exchange the second and third figure, but in the limit P 2 → P 1 the result does not change.
When a function f (z 1 , · · · , z n ) depends on n > 1 complex variables and there is a unique threshold, the singularities (solutions of 1/f = 0) are generically a subspace S ⊂ C n of codimension two and the branch subspaces V have codimension one, with S = ∂V. Thus, there are still two ways to analytically continue the function from C n \V beyond S to a neighborhood A of V. Again, the average continuation is half the sum of the two.
In the presence of several thresholds, we have several subspaces V. Their intersections give new regions A. To reach the intersection of two subspaces V we must perform two average continuations in different variables. It is easy to check that the result is independent of the order of the continuations. For example, let n = 2 and V i = {(z 1 , z 2 ) : Re[z i ] > 0}, i = 1, 2. Then we can reach the intersection V 1 ∩ V 2 either by first average-continuing in z 1 and then in z 2 , or vice versa, but the result does not change. The argument easily extends to multiple intersections.
We define the average continuation recursively. Consider an arbitrary diagram G. Deform the masses, so that the LW thresholds (2.25) are all distinct. Let G i 1 denote the result of the average continuation in some analytic region A i 1 , already reached, with nonvanishing widths . In the zeroth step, we take the result G 0 of the loop integral in the main region A 0 . We want to reach a new analytic region A i 2 above some LW threshold P . Redefine the external momenta p 1 , · · · , p n so that P reads p 0 j = p 2 j +M 2 for some p j and some combinationM of masses. Assume that an open-ball neighborhood U P of P belongs to A i 1 , apart from the points of the half line O P with Re[p 0 j ] > 0, Im[p 0 j ] = 0 in the p 0 j complex plane, wherep 0 j ≡ p 0 j − p 2 j +M 2 . The average-continued function in Here and below, the dependence on p j and the other external momenta is understood. After the evaluation of G AV , the deformed masses are sent back to their original values and the result found in U P ∩ O P is extended to a neighborhood of O P by analytic continuation, which defines the analytic region A i 2 above the threshold, as explained before. The operations (4.4) must be applied to every LW threshold. The relation between the average continuation and the nonanalytic Wick rotation can be proved as follows. LetÃ i 2 denote the region identified by the condition D pinch = 0, where D pinch is given by formula (2.16). The behavior of the loop integral around the pinching singularity insideÃ i 2 is, after integrating on the loop energies by means of the residue theorem,
where D k,q is the integration domain on the loop space momenta k and q. The denominator is a complex function, so its vanishing amounts to two real conditions. Write where f j and g j are real functions. When D k,q is deformed, the regionÃ i 2 is deformed as well. We have to arrange the domain deformation so as to squeezeÃ i 2 onto O P . Note that the deformed integration domain may depend on the external momentum p, as discussed in the previous section. We denote it by D def k,q (p). Referring to figure 16, we arrange D def k,q (p) so thatÃ i 2 turns into half a stripÃ def i 2 of thickness 2σ centered in O P . The parameter σ will later tend to zero, to complete the domain deformation and squeezeÃ def i 2 onto O P . When the loop momenta span the domain D def k,q (p), the external momenta span the regionÃ def i 2 . We can use this map D def k,q (p) →Ã def i 2 to make a change of variables such that −Re[p 0 j ] + f j (k, q) = τ and Im[p 0 j ] + g j (k, q) = ση. Then, in spacetime dimensions D greater than or equal to three, the integral gets the form where ∆ > 0 and h is regular at τ = η = 0. We understand that the integral over the remaining variables has already been made. When σ tends to zero, we obtain where P denotes the principal value and sgn is the sign function. This is the result of the nonanalytic Wick rotation.
To perform the average continuation, we replace p 0 j by p 0 j + iδ, with δ real and small. Then, we first take σ to zero keeping |δ| > 0 (which amounts to squeezing the regionÃ def i 2 onto O P ). At a second stage, we study the limits δ → 0 + and δ → 0 − . So doing, we approach O P from above (Im[p 0 j ] > 0, δ → 0 + ) and from below (Im[p 0 j ] < 0, δ → 0 − ). Since |δ| is small, the integral (4.5) becomes Averaging the two outcomes, we get (4.6) again. Thus, the nonanalytic Wick rotation and the average continuation give the same results, as claimed.

JHEP02(2018)141
With multiple thresholds the conclusions are the same, as long as the threshold locations are distinct, as emphasized before. For two thresholds located in τ 1 and τ 2 , we have integrals of the form . (4.8) If τ 1 = τ 2 the distributions of the form δ(τ − τ 1 )δ(τ − τ 2 ) that would appear in the limits σ → 0, δ 1 → 0 ± , δ 2 → 0 ± , vanish. Note that they are multiplied by a power η 2 , which is not killed by the integral over η. The distributions are instead killed by the integral over η (or the averages over the limits δ 1 → 0 ± and δ 2 → 0 ± ), so in the end we remain with both in the case of the average continuation and in the case of the nonanalytic Wick rotation. The arguments and conclusions easily extend to D = 2 once the integrals over η that appear in formulas (4.5), (4.6), (4.7) and (4.8) are replaced by sums over the values η = −1 and η = 1.
We conclude this section by mentioning other integral representations of the average continuation, which will be useful for the proof of perturbative unitarity. For definiteness, we take a unique LW threshold P and assume that it is located on the real axis. We deform the integration domain D k,q to a D +def k,q such that the boundary curve γ of figure 3 is turned into a curve γ like the one of figure 17. Then we consider the loop integral obtained by replacing the domain D k,q with D +def k,q . Clearly, this integral representation allows us to move analytically from the portion of the real axis that is located below the intersection with γ to an interval I of the real axis above P , without encountering LW pinchings. Let J + denote the result of the loop integral calculated in I following this procedure. At a second stage, we make a mirror deformation D −def k,q , so as to obtain a picture where γ is turned into the reflection of γ with respect to the real axis. We calculate the loop integral in I and call the result J − . The integral representation of the average continuation in I is (J + + J − )/2. We can further deform the domains D ±def k,q so as to stretch I to the whole O P . The construction easily generalizes to LW thresholds that are not on the real axis and to multiple LW thresholds.

Average continuation in various dimensions
In this section we illustrate the average continuation in examples related to typical loop integrals. The first example is f (z) = ln z, with the branch cut on the negative real axis. The functions f ± (z) of the previous section are ln(z ± i ), so, by formula (4.1), the averagecontinued function turns out to be The imaginary axis divides the complex plane in two disjoint regions: the half plane Re[z] > 0 and the half plane Re[z] < 0. The superior function can be determined from the inferior function, but neither of the two can be analytically continued beyond z = 0. By comparison, the Feynman prescription gives ln(z − i ).
The difference continuation gives which may be written as iπθ(−z), where θ(z) is the complex θ function, equal to 1 for Re[z] > 0 and 0 for Re[z] < 0. Note that the function ln z with z → −p 2 is the value of the bubble diagram of a massless scalar field in four dimensions, apart from an overall factor and an additional constant. The Feynman prescription leads to ln(−p 2 − i ), while the average continuation leads to f AV (−p 2 ) = (1/2) ln(p 2 ) 2 [19]. If we squeeze the half plane Re[z] < 0 onto the negative real axis, formula (5.2) encodes the discontinuity of the amplitude of the bubble diagram, i.e. the sum of the two cut diagrams associated with it, which is proportional to f d (−p 2 ) = iπθ(p 2 ).
As a second example, consider the function f (z) = √ z. We find Here, the superior function cannot be determined from the inferior one, which vanishes. The inferior function can be trivially continued beyond z = 0, while the superior function obviously cannot. In three dimensions, the bubble diagram of a massless scalar does not give a logarithm, but 1/ −p 2 . The Feynman prescription leads to 1/ −p 2 − i . If we use property (v) of section 4 with g(z) = 1/z, f (z) = √ z and h(z) = 1/ √ z, we find f AV (z) = 0 for Re[z] < 0. Again, the difference continuation is proportional to the discontinuity of the bubble diagram.

Four dimensions
In the massive case, the bubble diagram of the standard scalar field in four dimensions leads to the well-known expression after renormalizing the divergent part. This function has branch cuts in p 2 = (2m) 2 . Switching to the dimensionless variable z = p 2 /m 2 , we are lead to study the function In figure 18 we show the plot of this function for z real, together with the plot of the difference continuation. The first plot has the typical form of the LW amplitudes around the LW pinching [7]. Basically, the average continuation turns the ordinary scalar field into a massive fakeon (see the next section for details), i.e. the massive version of the fake degree of freedom of ref. [19]. Now, consider the LW propagator The bubble diagram built with it has the LW threshold p 2 = 2M 2 . Again, in the Euclidean region |Re[p 0 ]| < |p| we can evaluate the loop integral straightforwardly by means of the Feynman parameters. We have the sum of four contributions

JHEP02(2018)141
where a, b = + or −. The functions r ++ and r −− can be analytically continued to the whole real axis, because they are not interested by LW pinchings. Renormalizing away the divergent part, their sum is equal to −ig(p 2 /M 2 )/(8π) 2 , where For p 2 < 0, the sum of r +− and The function f 0 (t) does not give the correct result for t > 0. Indeed, it is symmetric under t → −t and not analytic in t = 0 (which is not a LW threshold). We have to analytically continue f 0 (t) from t < 0 up to the LW threshold t = 2. Then we have to average-continue it beyond the LW threshold.
Observe that K(x, t) has four zeros in x, which are We have to concentrate on the interval 0 < t < 2. We see that Im[v(t)] does not vanish, while Im[u(t)] vanishes for t = 0 and only there. In that point, u is equal to 1/2, which belongs to the integration path 0 < x < 1. When t grows and crosses the value 0, two zeros, u(t) and u * (t), cross the integration path, while the other two remain far away. It is simple to analytically continue the derivative f 0 (t) beyond t = 0, because its integrand is meromorphic. We just have to add the residues of the poles that cross the integration path, which are equal to −2πiu (t) and 2πiu * (t). When we go back to the primitive, we obtain, on the real axis, the function which is indeed analytic for t < 2. At this point, it is easy to perform the average continuation above the LW threshold t = 2. Observe that the average continuations of f 0 (t) and 1/t are trivial, while the average continuation of the square root is zero, by formula (5.3). Thus, above the LW threshold we just have to drop the square root. The final result is (on the real axis) Its plot is shown in figure 19 and is very similar to the one of the massive fakeon shown in the left picture of figure 18. We can call it Lee-Wick fakeon.
Repeating the arguments for the more general LW propagator and focusing on r +− = r −+ (r ++ and r −− still being analytic on the real axis) we get where r = µ 2 /M 2 and σ = 2 √ r 2 + 1 + 2r. Finally, we study the nonanalytic Wick rotation of the Euclidean theory and compare it to the average continuation. We work with the propagator (5.5). The average continuation of the amplitude on the real axis is where the combinatorial factor 1/2 is included.
If we want to evaluate the amplitude by means of the nonanalytic Wick rotation, we have to make the calculation inside the regionÃ P of figure 3 and deform the integration domain on the loop space momentum as explained in section 3, tillÃ P squeezes onto O P , which is the portion of the real axis from P to +∞. The procedure is involved, but there are situations where the regionÃ P is sufficiently thin to make the actual deformation unnecessary. One such case is when the LW scale M is small. It does not even need to be so small, since in most formulas it is raised to the forth power.
A measure of the violations of analyticity and Lorentz invariance, which occur before the domain deformation, is given by the "distance" between the point P and the point P of figure 3, i.e. the difference between the values of p 2 in such two points. Expanding the difference for M small, we find The first term is Lorentz invariant, so it controls the violation of analyticity. The second term controls the Lorentz violation. We see that the Lorentz violation is much smaller than the violation of analyticity. Numerically, we should see an evident Lorentz violation  figure 20, which confirms what we have just said. From left to right, the three plots are |p| = 3, 2 and 1. In the first picture, where M = 1, the plots superpose below the minimum P , but evidently deviate from one another above P and P . In the second picture, which has M = 1/20, the agreement is good everywhere. In figure 21 we include the prediction of the average continuation for M = 1/20, which is the top graph. As predicted by the first term on the right-hand side of formula (5.11), we see a discrepancy in the interval 0 < p 2 2M 2 ∼ .005 (caused by the missing domain deformation) and agreement everywhere else.

Three dimensions
In three dimensions the bubble diagram built with the propagator (5.5) gives the functions Here it is more tricky to work with the integrands, so it is better to eliminate the Feynman parameters by evaluating the integrals explicitly. In the Euclidean region t < 0 we find It is important to take such functions exactly as they are written, because manipulations that look innocuous may actually conflict with the determinations of the square roots and the logarithms. We have chosen to write the formulas so that they have the correct expansions for t ∼ −∞.
By formula (4.3), the average continuation on the real axis is just the real part, which gives the bubble amplitude Re[f 0 (p 2 /M 2 ) + g(p 2 /M 2 )]. (5.12) As in four dimensions, the nonanalytic Wick rotation exhibits, before the domain deformation, violations of analyticity and Lorentz invariance. They are apparent at M = 1, and, say, |p| = 1, 1/2, 1/3, as confirmed by the left picture of figure 22. By the estimate (5.11), we expect that Lorentz invariance is quickly recovered at, say, M = 1/20, which is confirmed by the right picture of figure 22, where |p| = 1, 2, 3. Zooming in, it is possible to observe a slight discrepancy around p 2 = 0, which is the violation of analyticity due to missing domain deformation and estimated by the first term of (5.11).
By applying formula (5.12), we can compare the results for M = 1/20 with the ones of the average continuation. This gives figure 23. Again, we see small discrepancies between P and P , due to the missing domain deformation, but agreement below P , where no domain deformation is required, and above P , where the effects of the domain deformation are negligible.

Two dimensions
In two dimensions the bubble diagram with propagators (5.5) gives, in the Euclidean region t < 0, a result proportional to the sum f 0 (t) + g(t), where As before, the integrand of f 0 has four singularities on the imaginary axis of the complex x plane. Two of them cross the x integration path when t varies from negative to positive values, while the other two do not cross the integration path. Since the singularities are poles, the difference f (t) − f 0 (t) for 0 < t < 2 can be easily calculated by summing the two residues, multiplied by 2πi. We find .
Then, f AV (t) = f (t) on the whole real axis. Indeed, we know that the average continuation of the function ∼ 1/ √ t is zero below t = 0. From the point of view of the nonanalytic Wick rotation, the two-dimensional models are a bit different from the models in dimensions greater than or equal to three. The reason is that in two dimensions the LW pinching occurs only at the boundary of the regionÃ P of figure 3, but not inside. The result of a loop integral in O P is Lorentz invariant and analytic even before making the domain deformation. The only Lorentz violation we find in the intermediate steps is due to the fact thatÃ P extends to P . To recover Lorentz invariance, it is sufficient to ignore the function found insideÃ P below P and analytically extend the function found inÃ 0 from P to P .
We can show these facts numerically, by plotting the results of the calculations for real p 0 around the points P , P , with various values of |p|. In figure 24 we see four vertical lines. The first three, from the left to the right, correspond to |p| = 3, 2, 1, with M = 1. Their locations are those of the point P . We see that each pair of plots agree both below the smaller P and above the larger P .
The forth vertical line of figure 24 corresponds to the result of the average continuation. We see that the nonanalytic Wick rotation with no domain deformation and the average continuation agree both below P and above P , even if M is not small with respect to |p|. In conclusion, a great simplification occurs in two dimensions, where the domain deformation is not strictly required to make calculations by means of the nonanalytic Wick rotation. At the same time, we have learned how powerful the average continuation is, because it drastically reduces the calculational effort in all dimensions.

Fakeons
We have seen that the average continuation is a simple operation to overcome branch points. Then, it is natural to inquire what happens if we apply it to a physical degree of freedom. Consider for example, the bubble diagram of ordinary scalar fields, which can be formally obtained by letting M tend to infinity in formula (2.2). The propagator just has the circled poles of figure 1. After taking → 0, the bubble loop integral has two branch points on the real axis at p 2 = (m 1 + m 2 ) 2 . The branch cuts are the half lines p 2 (m 1 + m 2 ) 2 on the real axis. An different from zero gives the familiar Feynman prescription, which displaces the branch cuts a little bit from the real axis and thereby allows us to define the loop integrals above the thresholds by analytic continuation from the segment p 2 < (m 1 + m 2 ) 2 to the half lines p 2 (m 1 + m 2 ) 2 . The displacements in the bubble diagram and its conjugate diagram are symmetric with respect to the real axis. This originates the discontinuity of the amplitude and, ultimately, the propagating degree of freedom. After subtracting the ultraviolet divergence, the diagram gives, in the massless case m 1 = m 2 = 0, where we have included the combinatorial factor 1/2.
The average continuation can be viewed as an alternative prescription to define the loop integral above the thresholds. If we forget about , by setting it to zero from the start, we can still define the amplitude unambiguously above the thresholds by means of JHEP02(2018)141 formula (5.1), in which case the result becomes (for p real) The discontinuity is absent, so we have no propagating degree of freedom. Equivalently, we can say that we have a fakeon, a fake degree of freedom. The average continuation makes the physical field disappear from the spectrum. At the level of the Feynman rules, the fakeon prescription can be formulated as follows. We replace the propagator 1/(p 2 − m 2 ) with [19] which coincides with (5.9) apart from the notation, and let E tend to zero at the very end.
The limit E → 0 is regular, since it is just a prescription for the propagator. The results of this paper apply to the theories whose elementary fields have free propagators that contain: (i) ordinary poles, treated by means of the Feynman prescription (with infinitesimal widths ); (ii) LW poles, with finite LW scales M ; (iii) fakeons, defined by means of the prescription (6.3), with infinitesimal LW scales E.
The widths must tend to zero first and the LW scales E must tend to zero last. At E > 0 we have a LW model, because the poles of type (iii) are just like the LW poles of type (ii). In that case, we make the computations by means of the nonanalytic Wick rotation or the average continuation. The results of the next section ensure that the theory is perturbatively unitary for → 0 at E > 0. If we let E tend to zero at the very end, perturbative unitarity is preseved, since it holds for every nonzero E.
We can retrieve the fakeon (6.3) from the results of the previous section, by taking the limit M → 0. For example, if we let M tend to zero in formula (5.10), we get −i times (6.2), which is correct, since for M → 0 the propagator (5.5) is the usual scalar propagator 1/p 2 endowed with the fakeon prescription (6.3). In three dimensions we can take the limit M → 0 of formula (5.12), which gives θ(−p 2 )/(16 −p 2 ).
While the LW degrees of freedom (ii) require higher derivatives and have finite LW scales M , the fake degrees of freedom (iii) can be introduced even without higher derivatives and have infinitesimal LW scales E → 0. Yet, there is not a deep difference between the two. In this respect, recall that the numerators of the propagators, such as the one of (6.3), are not important in the study of the LW pinchings. From now on, we call fakeons both the LW degrees of freedom (ii) and the fake degrees of freedom (iii). We may speak of fakeon thresholds, instead of LW thresholds, fakeon scales, and so on. We call fakeon theories the theories that involve fakeons (of LW type or not) besides ordinary physical degrees of freedom. Every result of this paper applies to the most general fakeon theory in dimensions D greater than or equal to 2.

JHEP02(2018)141
Observe that if we plan to take M , or E, to zero, the nonanalytic Wick rotation simplifies enormously, because there is no need to make the domain deformation. A quick way to see this is provided by formula (5.11), which gives an estimate of the analyticity violations and the Lorentz violations that occur prior to the domain deformation. Clearly, they both disappear in the limit M → 0. A more detailed argument can be provided by means of formula (2.20). Assume that we may have a LW pinching, i.e. n = n + + n − > 0. The pinching condition D pinch = 0, which defines the regionsÃ i , i = 0, implies We see that the vertical sizes of the regionsÃ i , i = 0, are bounded by nM/ √ 2, which tends to zero for M → 0. This means that all the regionsÃ i , i = 0, squeeze onto the real axis in that limit. Thus, the fakeons with E, M → 0 do not need the domain deformation.

Perturbative unitarity
In this section we derive the cutting equations and prove that the fakeon theories are perturbatively unitary to all orders. We assume that the Lagrangian is local and Hermitian.
Writing the S matrix as S = 1 + iT , the unitarity relation SS † = 1, which is equivalent to T − T † = iT T † , can be expressed diagrammatically by means of the so-called cutting equations [9][10][11][12], which relate the discontinuities of the amplitudes to sums of "cut diagrams". The cut diagrams are built with the usual vertices and propagators, plus their Hermitian conjugates, as well as "cut propagators". The cut propagators play a crucial role, because they tell us which degrees of freedom are propagated by the theory. Precisely, they encode the key completeness relation, which allows us to derive the unitarity equation SS † = 1 from the cutting equations. If ghosts are present, the cutting equations are still meaningful, but lead to a pseudounitarity equation instead of SS † = 1.
We want to prove that the fakeon models admit a physical subspace V of states and unitary cutting equations. This means that, if we project the initial and final states |α , |β onto V , only states |n belonging to V propagate through the cuts of the cutting equations. In other words, the completeness relation Obviously, we cannot demand unitarity for arbitrary complex external momenta, because the physical momenta are real. Therefore, we derive cutting equations that hold in a neighborhood U R ⊂ P of the subspace of real momenta and conclude that, thanks to them, the S matrix is unitary for real (on shell) external momenta. Note that the cutting equations also hold off shell. Figure 25. ACE propagators.

JHEP02(2018)141
We can assume that the LW scales M are arbitrary and different from zero. Once perturbative unitarity is proved in that case, it also follows for evanescent LW scales E, as long as they tend to zero after the widths .
The strategy of the proof is as follows. We first derive more general versions of the cutting equations that hold when the external momenta belong to the Euclidean region and the widths are nonvanishing. Then, we extend the validity of those equations to U R ∩ A 0 by analytic continuation and prove that they have the expected, unitary form in the limit → 0. Third, we average-continue the generalized cutting equations of U R ∩ A 0 to U R ∩ A i , i = 0, at = 0. Finally, we show that, in the limit → 0, the equations have the correct unitary form in every U R ∩ A i .
We begin by recalling an important tool that we use in the proof, i.e. the algebraic cutting equations.

Algebraic cutting equations
The algebraic cutting equations [12] are particular polynomial identities associated with Feynman diagrams. Let . . N , denote N sets made of four variables each. An abstract marking, called polarity and specified by the superscripts + or −, is assigned to these variables. We say that σ + i , τ + i (resp. σ − i , τ − i ) are positive (negative) polar numbers and use them to define the propagators Consider a Feynman diagram G with I internal legs and V vertices. We may assume that G is connected. Equip the G internal legs with orientations. We say that a curve is oriented if the orientations of all its legs are coherent. We say that a loop, i.e. a closed curve, is minimal if it is not the union of two loops that have a vertex in common.
Assign an independent energy to each internal leg and assume that it flows according to the leg orientation. Then, impose the energy conservation at each vertex, with zero energies on the external legs. This leaves L = I − V + 1 independent energies e 1 , . . . e L . We can arrange the orientations and the energies so that the flow of each e i defines an oriented minimal loop and the energy flowing in each internal leg is a linear combination of e 1 , . . . e L with coefficients 0 or 1. In this case, the diagram is said to be oriented.
Build variants G M of an oriented diagram G by marking any number of vertices. We define the value P M of G M by means of the following rules. Give value 1 to each unmarked vertex and value −1 to each marked vertex. Assign the propagators shown in figure 25 to the internal legs of G M , where the dots denote the marked vertices. Then, P M is the product of the values associated with the vertices and the propagators.
The algebraic cutting equation associated with G is the polynomial identity where P G is a linear combination of polarized monomials. A polarized monomial is a product of polar numbers, one for each internal leg, where at least one loop γ is polarized. We say that γ is polarized if the polar numbers associated with the legs of γ are arranged so that, moving along γ, the polarization flips if and only if the leg orientation flips.
The main virtue of the identity (7.4) is that it isolates the terms (those collected on the right-hand side) that do not contribute to the diagrammatic cutting equations. Indeed, in typical applications the polarity of a polar number refers to the position of its poles with respect to the integration path on the loop energy. A polarized loop is a product of polar numbers whose poles are all located on the same side. Letting tadpoles and nontrivial numerators aside, which can be treated with little additional effort [12], if we apply the residue theorem to perform the integral on the energy of a polarized loop, the result is zero.
To give a few examples, consider the diagrams of figure 26. The oriented loops of the third diagram are 123 and 34. Instead, 124 is a nonoriented loop. Equipped with polar monomials such as σ + 1 σ + 2 τ + 3 , σ − 3 τ − 4 and τ + 1 σ + 2 σ − 4 , respectively, these loops become polarized. Examples of polarized monomials for the third diagram are σ + 1 σ + 2 τ + 3 σ − 4 and σ + 1 σ + 2 σ − 3 τ − 4 . The polynomial identities (7.4) associated with the diagrams of figure 26 are where the polarized monomials on the right-hand sides have been replaced by zeros, since in the end they do not contribute to the diagrammatic cutting equations.
The algebraic cutting equations are more general than the usual diagrammatic cutting equations that are met in physics, in the sense that no particular assumptions are made about the polar numbers, apart from their polarity assignments. In the usual applications to quantum field theory, z i are the ordinary propagators and w i are their complex conjugates. Moreover, u i and v i are the cut propagators, i.e. distributions of compact support, typically theta functions that multiply delta functions. Here it is not necessarily so. For example, we are free to keep the infinitesimal widths of the Feynman prescription different from zero and arbitrary. Being able to work at = 0 is crucial to prove the perturbative unitarity of the fakeon models.  integration path gets pinched. We call this kind of pinching LW selfpinching, since it does not involve different propagators, but the poles of the same (cut) propagator. A mirror picture with respect to the imaginary axis is obtained for v. Observe that when tends to zero the standard poles of the cut propagators also pinch the integration path. We call that pinching standard selfpinching.

JHEP02(2018)141
To describe the LW selfpinching more clearly, it is convenient to start from different polar numbers, located in more usual positions, as shown in figure 29. Specifically, we take where Ω 1 , Ω 2 have negative imaginary parts, together with τ ± = −(σ ∓ ) * . Now the polar number σ + has three poles located in the fourth quadrant, while the polar number σ − has three poles located in the second quadrant. For example, making the M dependence explicit by writing Ω ± (p, M ) = p 2 + M 2 ± , we can set Ω 1 (p) = Ω − (p, M ) and Ω 2 (p) = Ω − (p, M ) for some real M = M . For the arguments that follow, it may also be convenient to pick a different M for every propagator.
If we keep the definitions (7.6) and take the real axis as the integration path for the energies, we can derive the algebraic cutting equations (7.4) using the polar numbers (7.9), by applying the Feynman rules of the previous subsection. When we integrate on the loop momenta, the right-hand side drops out, which leads to the diagrammatic cutting equations where the sum is over the properly marked diagrams G M , i.e. the diagrams that contain at least one marked vertex and one unmarked vertex. The diagramḠ is the one with all marked vertices. Figure 29. Poles of the half propagators (7.9).

JHEP02(2018)141
We have taken nonderivative vertices, so far, but the arguments also work when the vertices are polynomials of the momenta and the free propagators have nontrivial polynomial numerators. We stress once again that the equations (7.10) that we obtain are more general than the usual cutting equations, since the widths do not need to be small or tend to zero, but are completely arbitrary.
As long as the polar numbers are (7.9), the Wick rotation is straightforward. Then, however, the cut propagators do not simplify as in (7.7) and do not reduce to the expected form (7.8) when → 0. We must migrate Ω 1 (p) to Ω + (p, M ), which is equivalent to complexify M and deform M 2 continuously into −M 2 . During the migration, Ω 1 crosses the real axis. To keep the algebraic cutting equations valid, Ω 2 cannot cross the integration path, since the definition of polarity refers to the positions of the poles in p 0 with respect to it. Thus, we have to deform the integration path so as to avoid the crossing. This operation, applied on σ + , leads to the first picture of figure 27. It leads to the second picture of figure 27 when it is applied to τ − = −(σ + ) * . When we apply it to the cut propagators u and v, we must take into account that the LW pair of σ + and the LW pair of (σ + ) * remain on opposite sides of the integration path, which leads to figure 28 and its reflection with respect to the imaginary axis. This is the reason why we cannot drop the LW pairs from the difference u = σ + − (σ + ) * so quickly. First, we have to make one LW pair cross the integration path. Once it is on the other side, it does "annihilate" the other pair. However, the crossing leaves a remnant (the contributions of a pair of residues), which must be taken into account. To prove perturbative unitarity we need to show that such a remnant does not contribute to the cutting equations.
Observe that the crossing only concerns the cut propagators. In uncut propagators, the migration of the poles Ω 1 (p) just returns the right result, shown in figure 1, and no selfpinching occurs. For this reason, the left-hand side of the cutting equation (7.10) goes directly to its correct, final form. Only the right-hand side needs a detailed analysis.
Consider a properly marked diagram G M . Assume that the cut propagators are n + 1 and depend on n loop momenta (the most general case being a straightforward generalization of this one). Each cut propagator u = σ + + τ − and v = σ − + τ + receives contributions from LW selfpinchings and standard selfpinchings. We decompose G M as a sum of terms where each cut propagator involves either of the two. We analyze such terms one by one, starting from the terms that involve only LW selfpinchings.

JHEP02(2018)141
Integrate on the n loop energies k 0 i by means of the residue theorem and take M 2 → −M 2 in n cut propagators. This operations give n conditions of the form k 0 i =ω i (k i ), which eliminate the loop energies k 0 i . At this point, the contribution of the LW selfpinching due the last cut propagator has the form where D ± pinch are deformed versions of the denominators D pinch of equation (2.16). The deformations depend on M and are such that D ± pinch → D pinch when M 2 → −M 2 . Moreover, they make D ± pinch vanish on opposite sides of the integration path. After the integrations on k 0 i , the integration path has actually disappeared, so formula (7.11) can be read as it stands. When we finalize the migration of Ω − (p, M ) into Ω + (p, M ) by taking the limit M 2 → −M 2 , the difference (7.11) gives zero, because we are working in the Euclidean region, where the loop space momenta are integrated on their natural real domains and the condition D pinch = 0 has no solutions. We recall that, indeed, D pinch = 0 is the condition for having a LW pinching, which defines the regionsÃ i , i = 0. Now, consider the terms where only standard selfpinchings occur. Those are the expected terms, the only ones that should survive at the very end. Indeed, the differences (7.11) give (7.7) in this case.
Finally, consider the mixed selfpinching, i.e. the terms where the contributions of some cut propagators come from LW selfpinchings and those of other cut propagators come from standard selfpinchings. Recall that the LW selfpinching occurs when we complete the migration of Ω − (p, M ) into Ω + (p, M ) by taking M 2 → −M 2 . Instead, the standard selfpinching occurs when we take → 0. If we are willing to let the widths disappear at the end, the argument used for the terms with only LW selfpinchings can be applied with straightforward modifications and leads to the conclusion that the contributions of the mixed selfpinchings vanish in the limit → 0. For various arguments that follow, however, it is necessary to keep = 0. There, we have generalized cutting equations that contain extra contributions, which must be taken into account for the extension of the proof beyond the Euclidean region. For example, consider the case where the contributions of the first n cut propagators come from LW selfpinchings and those of the last cut propagator come from a standard selfpinching with mass m. We integrate on the n loop energies as before and complete the migrations M 2 → −M 2 . At the last step, we obtain an integrand proportional to an expression of the form (7.11), where the denominators D ± pinch are equal to (2.20) with r = 1 and imaginary parts ∓i attached to the squared mass m 2 . Clearly, (7.11) does not vanish in this case until we take → 0.
Summarizing, the expected, unitary cutting equations hold in the Euclidean region for → 0. The cut propagators can be effectively replaced by (7.8) in that limit and the LW degrees of freedom do not propagate through the cuts. Moreover, generalized cutting equations hold at = 0.

Perturbative unitarity in the other regions
The next step is to extend the validity of the generalized cutting equations by analytic continuation from the Euclidean region to the intersection U R ∩ A 0 . Then we have to reach the other regions U R ∩ A i , i = 0, by means of the average continuation. In both cases, we must prove that the generalized cutting equations reduce to the expected, unitary cutting equations in the limit → 0. We assume that the masses are arranged so that the LW thresholds are all distinct.
We have seen that the generalized cutting equations in the Euclidean region have corrections C(p, ) for = 0, due the mixed selfpinchings, where p are the incoming external momenta. The reason why they vanish for → 0 is that D pinch never vanishes in the Euclidean region.
The first extension away from the Euclidean region is straightforward. At = 0 the standard branch points are displaced from the real axis. Moreoveor, we know that we can deform the integration domain on the loop space momenta so as to avoid the LW pinchings everywhere in A 0 . Once we do that, we can analytically continue the generalized cutting equation (7.10) from the Euclidean region to U R ∩ A 0 by keeping = 0 and moving along the real axis. Then, the corrections C(p, ) still vanish when we take the limit → 0, because D pinch never vanishes.
When we attempt to analytically continue the cutting equation (7.10) above an LW threshold P , we find that it cannot be done in a unique way. Averaging the two independent ways of doing it, we can prove perturbative unitarity in the regions U R ∩ A i , i = 0.
Specifically, we make the two domain deformations D k,q → D +def k,q and D k,q → D −def k,q explained at the end of section 4. Applying the deformations on the entire cutting equation (7.10), we obtain two deformed versions of it.
In the case of the deformation D k,q → D +def k,q , we denote the deformed versions of the diagrams G,Ḡ and G M by J + ,J + and J +M , respectively. In the case of the deformation D k,q → D −def k,q , we denote them by J − ,J − and J −M . In each case, we obtain an integral representation of the cutting equation (7.10) in some interval I of the real axis above P and we can reach I by analytic continuation from the Euclidean region without encountering LW pinchings. Since D pinch never vanishes in I, the corrections C(p, ) still vanish for → 0.
Note that the left-hand sides J ± +J ± of the deformed cutting equations are no longer real, because the integral representations of J ± andJ ± have the same (complex) deformed domains D ±def k,q . By construction we haveJ ± = (J ∓ ) * . When we average the two deformed cutting equations, we obtain the cutting equation that holds above the LW threshold. The average of the left-hand sides gives where (J + + J − )/2 is the average continuation of G and (J + +J − )/2 is the average continuation ofḠ. The average of the right-hand sides has the expected form for → 0, since the contributions C(p, ) drop out in that limit. The conclusion holds in the neighborhood of every I ⊂ U R ∩ A i , so it also holds in the whole U R ∩ A i . Applying this procedure to each LW threshold at a time, we reach JHEP02(2018)141 every U R ∩ A i , i = 0. When anomalous thresholds are met, there are multiple ways to circumvent them, which correspond to multiple options for the deformations, as described at the end of section 3. Each option can be used to average-continue the cutting equations as described above. The corrections C(p, ) vanish for → 0 in every case.
In the end, the cutting equations have the expected unitary form in all the regions A i for → 0. This concludes the proof that the fakeon models are perturbatively unitary to all orders. Note that it would be much more difficult to make the extension to A i , i = 0, using the nonanalytic Wick rotation. This shows once more the power of the average continuation, a very simple operation that allows us to make a number of manipulations that otherwise would be very cumbersome.

Remarks
Before concluding this section, we comment on the resummation of the perturbative series and its effects on the unitarity equation SS † = 1. We recall that the LW poles of the free propagators (2.1) are located symmetrically with respect to the real axis. This is important for the proof of perturbative unitarity, because the contributions of complex conjugate LW poles compensate each other. However, the exact two-point functions may lose the symmetry just mentioned, because the resummations may give widths to the standard poles and the LW poles, and change their masses. This is no source of concern, because that symmetry, which is helpful to see unitarity at the perturbative level, plays no role after the resummations.
Once we have derived the diagrammatic cutting equations (7.10) and projected the external states onto V , we have the completeness relation (7.1) and the unitarity equations (7.2). At a first stage, let us ignore the resummations that affect the standard poles and concentrate on those that affect the LW poles. Then the states of V stay the same and the unitarity equations (7.2) remain valid. These types of resummations just act internally to the correlation functions associated with α|T |β , α|T † |β , α|T |n and n|T † |β . At a second stage, we perform the resummations that affect the standard poles. Some physical particles may acquire widths and decay, and so disappear from the physical spectrum at very large distances. Since they still propagate through the cuts of the cutting equations, the S matrix is no longer unitary in a strict sense, although it remains perturbatively unitary.
In other words, when we resum the perturbative expansion, the LW sector does not affect unitarity. Yet, some physical poles may get nonvanishing widths, pretty much like the muon in the standard model. In this respect, the fakeon models behave as an ordinary model.
If the Lagrangian is Hermitean, the results of the next section ensure that its renormalization is also Hermitean, so the denominators of the renormalized propagators obtained by including the counterterms still have the structure displayed in formula (2.1), with pairs of complex conjugate poles, besides the physical poles.

JHEP02(2018)141 8 Renormalizability
Commonly, higher-derivative theories are thought to have an enhanced power counting, because the propagators fall off more rapidly at high energies. However, the usual rules of power counting just work in Euclidean space, while in Minkowski spacetime it is much more difficult to have control on the ultraviolet behaviors of the Feynman diagrams. Everything is fine if the Minkowski formulation of the theory is analytically equivalent to the Wick rotated Euclidean one, which happens for example when the free propagators just have poles on the real axis. A fakeon model does not have this property, to the extent that the Minkowski version is plagued by nonlocal, non-Hermitian counterterms [30]. At the same time, we know that the Wick rotation of the Euclidean version of a fakeon model is not analytic everywhere, so we have reasons to worry that the nice renormalizability properties of the Euclidean version may not be fully inherited by the nonanalytically Wick rotated theory.
In this section we overcome these worries, by proving that the renormalization of a fakeon model is still local and actually coincides with the one of its Euclidean version. We give two arguments, the first one based on the average continuation and the second one based on the nonanalytic Wick rotation.
The first argument is straightforward. Once we have subtracted the divergences of the Euclidean theory, the amplitudes are convergent in the Euclidean region. We know that we can unambiguously reach every other region from there. The analytic continuation of a convergent function is obviously convergent. The same holds for the average continuation, which is made of two analytic continuations. This implies that the amplitudes are fully convergent in every analytic region A i .
The second argument requires a bit more work. The rules of power counting of the Euclidean theory trivially extend from the Euclidean region to the main region A 0 , since the Wick rotation is analytic there. So, we just need to concentrate on the other regions A i , i = 0. Let us start from the regionsÃ i , i = 0, which are defined as the solutions of the conditions D pinch = 0 with real loop space momenta, D pinch being given by (2.20). As we know, the relative sign in front of the frequencies of (2.20) is necessarily positive, otherwise no pinching occurs. Assume that the external momenta p belong to a compact connected open subset S p ⊂ P that contains an open subset of the Euclidean region. Formula (2.20) makes it clear that the condition D pinch = 0 cannot be satisfied in S p for arbitrarily large |k i | and |q j |. Thus, the solution identifies a compact subset C k,q of the domain D k,q of the loop space momenta.
Recall that the loop energies k 0 i are gone after applying the residue theorem. Split the integral on D k,q as the sum of the integral on a compact subset C k,q ⊃ C k,q plus the integral on D k,q \C k,q . Clearly, the integral on C k,q is not interested by ultraviolet divergences. On the other hand, the integral on D k,q \C k,q may be ultraviolet divergent, but it is not interested by the LW pinching. This means that it admits an analytic Wick rotation, which makes its ultraviolet divergences equal to those of its Euclidean version. Observe that the Euclidean loop integral is reachable analytically while remaining inside S p , since S p is chosen to contain an open subset of the Euclidean region. Thus, once the JHEP02(2018)141 (Euclidean) divergences and subdivergences are subtracted, the loop integral is convergent in S p . Since S p is arbitrary, the subtracted integral is convergent everywhere in P.
So far, the integration domain D k,q is still undeformed, because we have been working in the regionsÃ i . Now we have to perform the domain deformation to go from the regions A i to the regions A i . We can make it so that the deformed C k,q remains always compact. Applying the argument above to every deformed D k,q , we see that the final result is convergent in every region A i .
We conclude that the nonanalyticity of the Wick rotation does not conflict with the renormalization of the fakeon models, which coincides with the renormalization of their Euclidean versions. In particular, the locality of counterterms and the usual rules of power counting hold. This proves that the fakeon models that are renormalizable do reconcile unitarity and renormalizability.

Conclusions
In this paper we have studied the fakeon models, which contain ordinary physical particles and fakeons, i.e. fake degrees of freedom. An important subclass are the Lee-Wick models, which have higher derivatives. Fakeons can also be introduced without higher derivatives, by means of a suitable quantization prescription.
Formulating the models by nonanalytically Wick rotating their Euclidean versions, we have shown that they are consistent to all orders. In particular, we have studied the LW pinching and the domain deformation in arbitrary diagrams.
The S matrix of the fakeon models is regionwise analytic. Different analytic regions A i are related by the average continuation, a powerful operation that allows us to simplify numerous derivations. The average continuations of various functions that are frequently met in four, three and two dimensions have been computed and compared numerically to the results of the nonanalytic Wick rotation, confirming that the two operations give the same result.
We have proved that the fakeon models are perturbatively unitary to all orders. The strategy of the proof was to first use the algebraic cutting equations to derive generalized versions of the diagrammatic cutting equations that hold in the Euclidean region at = 0. Then we have shown that the equations can be analytically continued to the main analytic region A 0 and average-continued to the other analytic regions A i , i = 0. Finally, we have proved that they reduce to the expected, unitary cutting equations when the widths tend to zero.
Another good property of the fakeon models is that they have the same renormalization as their Euclidean versions have. This makes them viable candidates to explain quantum gravity. We recall that while the LW models of quantum gravity [19,28,29] are superrenormalizable, the fakeon models of quantum gravity can be strictly renormalizable [19]. At present, the best candidate to explain quantum gravity is a fakeon theory in four dimensions whose Lagrangian density contains the Hilbert-Einstein term R, the cosmological term and the terms R µν R µν , R 2 [19]. It is the unique model whose gauge coupling is dimensionless. It has all the features we expect apart from one: a nonvanishing JHEP02(2018)141 cosmological constant, which may predict a small unitarity anomaly in the universe. The classical action of this theory coincides with the one considered in refs. [15][16][17][18] and more recently refs. [36,37], but its quantization and physical predictions are completely different, because the would-be ghosts have been replaced by the fakeons. Strictly unitary superrenormalizable models can also be built [19], but their features makes them less realistic. In the end, the fakeon models have all the features that we require to include them into the set of the physically acceptable theories.