Unitary evolution of a pair of Unruh-DeWitt detectors calculated efficiently…

Unitary evolution of a pair of Unruh-DeWitt detectors calculated efficiently to an arbitrary perturbative order

Kamil Brádler Department of Mathematics and Statistics, University of Ottawa, Ottawa, Canada

Unruh-DeWitt Hamiltonian couples a scalar field with a two-level atom serving as a particle detector model. Two such detectors held by different observers following general trajectories can be used to study entanglement behavior in quantum field theory. Lacking other methods, the unitary evolution must be studied perturbatively which is considerably time-consuming even to a low perturbative order. Here we completely solve the problem and present a simple algorithm for a perturbative calculation based on a solution of a system of linear Diophantine equations. The algorithm runs polynomially with the perturbative order. This should be contrasted with the number of perturbative contributions of the scalar theory that is known to grow factorially.

Speaking of the model, a welcome collateral result is obtained to mechanically (almost mindlessly) calculate the interacting scalar theory without resorting to Feynman diagrams. We demonstrate it on a typical textbook example of two interacting fields for .

Key words and phrases:
Unruh-DeWitt detector, Isserlis’ and Wick’s theorem, interacting scalar model, Feynman diagrams

1. Introduction

An interaction model between a real scalar field and a two-level atom called Unruh-DeWitt (UDW) particle detector was conceived [1], improved [2] and studied in many physical situations [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. A UDW detector consists of a two-level atom linearly coupled to a scalar field. Research in this area typically focuses on a single UDW detector to probe the field or two UDW detectors to study bipartite correlations with a special interest in the vacuum entanglement. The calculation is of a perturbative character but to the author’s knowledge, the Unruh-DeWitt Hamiltonian has not been (dis)proved to be exactly solvable [15]. The main result of this paper is the next best thing. We also approach the problem perturbatively and provide an expansion to an arbitrary order of a coupling constant and so we can study the field-detector(s) interaction in the strong coupling regime. Crucially, unlike a typical situation in quantum field theory, we show that the number of perturbative terms grows polynomially with the perturbative order. This is in a sharp contrast with the textbook problem of the interacting scalar theory (and other more realistic models) where the number of terms grows factorially [16, 17, 18, 19]. Additionally, the perturbative series are typically asymptotic and this limits even the perturbative approach itself. As one reaches the perturbative order of approximately the inverse of the coupling strength the series summands start to increase. The series diverges [20, 21] and resummation techniques must be used [22]. In the case of two UDW detectors we show that there is no factorial explosion and the polynomial growth indicates no problems with convergence. Could it be then that the actual series for two UDW detectors we obtained here is summable or is the issue of divergence more subtle and persistent? This is an interesting open question probably waiting to be rigorously settled but the results presented in this paper will prove to be useful irrespective of the answer.

Surprisingly, a part of our UDW calculation turns out to be useful for a perturbative calculation of the scattering amplitudes in the interacting scalar theory. A typical approach for the high-order loop corrections is to enumerate all participating Feynman diagrams [23, 19]. Here we circumvent this step and instead use the statistics cousin of Wick’s theorem [24] known as Isserlis’ theorem [25] to directly calculate the perturbative contributions. This is possible due to an easy way of finding all (non-negative) solutions of a system of linear Diophantine equations Isserlis’ theorem leads to. As a result, the scattering amplitudes for the interacting scalar model can be found without the need to infer Feynman’s rules [26]. Feynman’s diagrams can, however, be easily recovered and we demonstrate the whole process on the textbook examples of the second order expansion of the and  models (for two interacting fields). This obviously does not remove or circumvent the fast growth of the perturbative contributions in these models but it seems that the approach coming from this direction is novel and makes the perturbative calculations (even to high orders) extremely streamlined. This may be surprising due to the fact that Diophantine equations have a well-earned reputation of being hard to solve.

In Section 2 we recall the definition of the Unruh-DeWitt detector, take two copies of it and introduce the problem we would like to solve. The solution is presented in Section 3 and is divided in three phases. In Section 4 we show how to apply our approach to perturbatively calculate the scattering amplitudes for the scalar model.

2. Unruh-DeWitt detector

A realistic UDW detector is described by the following term




is a real scalar field where is a four-momentum, is a detector’s window function ( is its proper time and are Minkowski coordinates), a smearing function, are the detector ladder operators, stands for the energy gap and for a coupling constant. The evolution operator for observer reads


where is the time-ordering operator, and similarly for observer .

Out task is to efficiently calculate the following unitary


for any , no matter how large, and for any input atomic state. We denoted


and so holds. Note that we set to be a delta function for convenience. The result of this paper is oblivious to the details of smearing, window functions or trajectories. The important expression is Eq. (4) and the fact that are bosonic in nature. Our final goal is to express -point correlators in terms of a polynomial number of 2-point correlators and not their calculation where these details are relevant.

The matrix elements of (4) (in the canonical basis) will be shorthanded as


where and a tensor product implied. So for the detector density matrix we get


where denotes Minkowski vacuum and is the detectors’ initial state.

3. Two Unruh-DeWitt detectors to an arbitrary perturbative order

This section contains our main result. We will derive for and show a trivial modification enabling us to calculate Eq. (7) for any canonical basis state . This, on the other hand, will lead to a ‘standalone’ unitary matrix that can be used to find for any detector input state and thus completely solving the problem. In the first step (Phase I) we tame the exponential number (in ) of summands of the core expression Eq. (7) (or (4)). By plugging the result into Eq. 8 we get again a polynomial amount of -point Green’s functions. We will be able to split them into a product of two-point correlation functions by constructively solving a system of linear Diophantine equations with a polynomial number of non-negative solutions in (Phase II) and calculate their multiplicities (Phase III).

Phase I

Let us recall the elements of Pauli operators algebra representing the two-level detector(s). They satisfy and also together with . That is, the atomic operators for the two detectors commute in the Unruh-DeWitt model.

Proposition 1.

For we obtain from Eq. (4)


and for we get


The initial state is a ground state: . There are two possibilities for the output states: and but taking into account the nilpotency property of we can also enumerate all possible ways the output states can be obtained. It is simply


for all . All other operator sequences result in zero. In our case we have two sets of ladder operators – for the and Hilbert spaces and so


In the first two lines is even and in the last two lines is odd. The sums should be understood as counting all possible strings of the sigma operators leading to a given output state on the left and we will now find the explicit expressions. Let’s call them nonzero strings. Given the desired output state, the main observation is that there are always only two sigma operators that do not destroy the previous nonzero string. This is because every nonzero string acting on can land in only one of the four possible states. For instance, if the state is the next sigma operator can be or to get another nonzero string. Similarly for the other three states and so we can represent all nonzero strings as a complete binary tree of the depth where we adopt the following convention: The root node is the initial state where the left offsprings correspond to the action of and the right offsprings correspond to . The nodes are labeled by the resulting state , see Fig. 1. In all the branchings that follow the left offsprings correspond to the action of and the right offsprings correspond to . Given the offspring node, the choice of either or in the branching is unambiguous.

Figure 1. All paths for the first two branchings starting from are depicted. The arrows indicate the action of the ladder operators . The resulting state is a new node.

Figure 2. A schematic depiction of the action of Eqs. (1) and (1) for branchings. As in Fig. 1 the arrows symbolize the action of the detector ladder operators . The left-pointing arrows correspond to the detector while the right pointing arrows are for the detector. The initial state is and the symbols and are the boson operators ‘left behind’ after the action of . So as we follow any path in the tree (up or down from ) we collect the boson operators thus forming an operator sequence in Eqs. (1) and (1). There are many paths with the same operator products and their multiplicities are the binomial coefficients derived in Proposition 1.

But this immediately provides the desired counting because we have just proved a bijection between the number of branches of a perfect binary tree of the depth and the number of non-zero strings of the same length. Consider Eq. (13a). We set the total number of branchings and giving us . The number of paths with left branchings (one for and one for making it even since we have to end up in ) out of the total branchings is and so is the number of nonzero strings in Eq. (13a). Identically, for (13b) we again set and and obtain of nonzero strings ( left branchings and right branchings). For Eq. (13c) we set and . Hence is the total number of nonzero strings. Finally, for Eq. (13d) we have and we now set . In this way we obtain the same number of nonzero strings as in the previous case.

The proposition statement is obtained by inspecting the LHS of Eqs. (1) and (1) where the boson operators accompany the corresponding ladder operators and the RHS consist of the products of the boson operators from nonzero strings. ∎


Fig. 2 illustrates the proof. It also illustrates the fact that the same analysis can be done for any initial state . Depending on the values of and the arrows in the binary tree will represent different ladder operators leading to the same counting and different boson operators on the RHS of Eqs. (1) and (1). In this way we are able to find all the unitary matrix elements and therefore we can act on any input state (pure or mixed) from a space spanned by .

Example 1.

Setting in Eq. (1) and in Eq. (1) (corresponding to the perturbative order ) we quickly reproduce the fourth order calculation in [14] (Eqs. (11)) and the result can be applied to a variety of situations [5, 6, 27, 8]. With the same ease we can obtain an output state for .

A non-trivial check is the trace of (8) being equal to one irrespective of the perturbative order. Case was verified in [14] and a straightforward calculation shows that it is true for as well111Note that the fourth order is an absolute must if one wants to calculate any entropic quantity to properly assess the importance of entanglement for quantum communication. If only the second order is calculated (almost always the case with a notable exception of [14]) then the component of (8) is zero and two of the eigenvalues are negative.. ∎

Phase II

We found the expressions and, crucially, the exponential growth in in Eq. (4) was managed such that the number of summands in Eqs. (1) and (1) increases polynomially with . However, upon inserting two such expressions into Eq. (8) we note that the exponential growth crept back. This is because of the following result:

Theorem 2 (Isserlis’ [25]).

Let be a Gaussian random variable satisfying . Then


where the sum goes over the products of bivariate expectation values .

Remark (notational).

In our case the Gaussian random variable will be a product of time-ordered free scalar fields and the expectation value will be taken wrt to Minkowski (i.e non-interacting) vacuum :


where and on the RHS is the minimalist notation we will be mostly using.


The theorem is also known as Wick’s theorem [24]. Wick’s theorem transforms time-ordered operator expression into a normal form [26] and Isserlis’ result is recovered upon taking a (free) vacuum expectation value. As a matter of fact, Isserlis’ theorem is more general since naturally there is no notion of time-ordering in Eq. (14) and so the theorem applies even for ‘ordinary’ products of free scalar fields. Of course, we are actually calculating the time ordered version as it appears in the definition of the sought after Green’s function.

The exponential growth lies in the presence of the double factorial in (14). For different fields in Eq. (15) there is not much to do but the situation is different for a constant number of fields. We first prove the following result whose application goes beyond the problem solved in this paper (see Sec. 4).

Theorem 3.

Let and . Then Green’s function can be written in terms of a polynomial number of products of two-point correlation functions:


where is the product multiplicity, and .

Definition 1.

Different products  will be called conformations and a conformation exponent.


At first sight there is nothing surprising about Eq. (16). Except for the double factorial the equation looks exactly like Eq. (14) and it must be that way. After all, it is its special case. The novelty lies in the polynomial number of summands for a constant number of fields (four in this case) and in a constructive (and simple) way of obtaining all the coefficients . The multiplicity factors are another crucial ingredient and they will be calculated in Theorem 6.


The solution of Eq. (16) can be found by solving the following system of four linear Diophantine equations for :


Linear Diophantine equations have zero or infinitely many solutions (counting both positive and negative ones). The numerical coefficients for the variables are so simple that we can easily guess a solution for (17a) (say and ) and so there are infinitely many solutions. It also means that there are finitely many non-negative solutions we are after. How do we obtain them in a systematic way? Since and for we fix (starting from ) and simply list all the admissible and . We continue by increasing by one and repeat the whole process and the same for . In the next step we move to Eqs. (17b)-(17d) by inserting all found triples and repeating the procedure until we find all admissible solutions in the second row. We insert the found solutions to the third and fourth row until we find all solutions. The result is a complete list of ten-tuples solving Eqs. (17).

How many different ’s are there? We don’t need an exact number (it is actually not that straightforward to count the solutions as we will see in Lemma 4) and a polynomial upper bound is enough222More precisely, we need to list all the solutions and once we have them we simply count them. What we don’t need for the purpose of this work is a closed expression giving the number of non-negative solutions of (17).. The number of solutions for Eq. (13a) depends on the parity of . For even and the number of nonnegative solutions is a triangle number and by incrementing by one up to we get the rest of solutions. Summing them up we obtain a total of


solutions. A similar analysis leads to the following number of solutions for an odd :


If all equations (17) were independent we would get a scaling for the number of solutions like , where . This is enough for the proof but it grossly overestimates the number of solutions. A more careful counting gives us the scaling. ∎


There is an extensive literature on how to solve a system of linear Diophantine equations [28, 29]. But the question is not really whether one can solve it but how fast (see, for example, [30]). It is likely that the simple form of our system does not even require any sophisticated method to solve it. One just systematically lists all the solutions as sketched in the proof.

Even a better estimate is given in the next lemma.

Lemma 4.

The number of non-negative solutions to (17) is upper bounded by


where .


To simplify the derivation, assume and if necessary let’s promote to be the closest (greater) even number. From all admissible it is yielding the greatest number of solutions so instead of (17) let’s start by studying the following simpler Diophantine system:


Eqs. (21) inherited an interesting structure from (17). Starting with (21b), every row shares exactly one variable (but always different) from all previous rows. From the last two lines of (21) we observe . Using the first two lines of (21) we also find and so we get


as the only consistent solutions. Then for any allowed pair of coefficients (allowed means for and ) there is one pair that can satisfy (22). This is because even though there are two other free variables, and , the permutation symmetry of (21) together with Eqs. (22) imply

Finally, looking at (21a), we find , where stands for the multiplicity of with the value , and so there is


solutions of (21). By upper-bounding all systems of linear Diophantine equations by (23) we get (20). This scales as . ∎

Example 2.

The number of conformations of is 17 and their exponents read



A calculation performed in [31] leads to a closed expression for the number of non-negative solutions of (17) for even to be


For we get from the previous example.

Phase III

Let’s find the expression for the conformation multiplicity factor . First an auxiliary lemma.

Lemma 5.

Let and be discrete sets where . Then there is


bijective functions where and where . For we set Eq. (26) to one.


If then (26) becomes which is known to be the number of bijections from a set to itself (the number of permutations). To make the notation more concise in following text we will use the definition of the falling factorial: . So Eq. (26) is written as .


The coefficient is simply a number of all possible domains whose cardinality is . Then for every domain there is codomains of the same cardinality. The coefficient comes from the number of permutations within each of codomains . ∎

Theorem 6.

For a given conformation exponent the multiplicity factor reads




We will use repeatedly Lemma 5 by identifying for from Theorem 3 in the following way:


So is the set of all ’s and is the set of all ’s. However, if we do more than one two-point Green’s functions we have to take into account that the sets and might have shrunk. This depends on whether in the preceding Green’s function we had or . The strategy we will follow here is to first count the multiplicity of as they are independent (meaning non-overlapping for different ’s) and then the multiplicities of and in this order(!).

The multiplicity of immediately follows from Isserlis’ theorem. Looking at Eq. (14) we see that if there will be identical products. Hence ()


as follows from Eq. (26)333Indeed, the double factorial follows either from Theorem 2 as stated or the counting in Lemma 5.. The factor of two in the binomial ‘denominator’ accounts for the two ’s in . Repeating this procedure for all we get in Eq. (28a). To get we then have to keep track of the set cardinality. To get (using the notation of Lemma 5) we set and notice that the cardinality of the discrete set of ‘ones’ is decreased by those contributing to . So . Similarly, the cardinality of the ‘twos’ is decreased by the number of elements contributing to . So and in (28b) follows. At this point it is evident how to proceed so let’s only derive which is next in line. We set and the set of ‘ones’ is now smaller by those contributing to but also to those from calculated in the previous step. Hence . The set of ‘threes’ was decreased only by the elements contributing to and so and Eq. (28c) follows. Going up to we arrive at the overall multiplicity factor, Eq. (27). ∎


We emphasize the importance of the order in which we follow the construction of . Other possibilities are certainly plausible but we believe the sketched one is the most intuitive.

Example 3.

(Important) Let’s find the multiplicities from Example 2. By setting and plugging into Eqs. (28) we get


where . Crucially (and this is a highly non-trivial check of both Theorems 3 and 6!), we get


According to Eq. (14) we have and so the RHS of (34) is the total number of products as it should be according to Isserlis’ theorem. ∎

Putting all pieces together

We use Proposition 1 to insert Eqs. (1) and (1) to Eq. (8) and get the output state components to the -th order. Let’s take a look at the component:


We identified and and in this form it is ready for Theorem 3 and onto Theorem 6. The number of terms in Eq. (3) is polynomial in and therefore in as well. By setting and taking the greatest , Lemma 4 gives us a tighter polynomial upper bound on the number of solutions compared to Theorem 3. We thus arrived at the main result of the paper.

4. Interacting fields in the scalar model to the second perturbative order

One of the simplest interacting theory is given by the Lagrangian


where is a real scalar field and a coupling constant. The model is not exactly solvable and its perturbative evolution, Feynman rules and much more (such as renormalization) are typically trained on  [26] or on a more physical model . The quantity of interest is again the -point Green’s function for interacting fields and its expectation value is calculated wrt to the interacting vacuum . Isserlis’s or Wick’s theorem cannot be applied for interacting fields but from the Lagrangian (or Feynman path integral) formulation [26] we know that the interacting Green’s function satisfies


Interestingly, for Theorems 3 and 6 can be used to an automatic calculation up to the second order for any as we will now illustrate for and .

Example 4 ( for ).

We expand the numerator of (37) as is customary [26, 32] and take the third term. It is proportional to


where on the RHS we use the notation of Theorem 3 by setting and . Similarly the second term of the denominator is proportional to


Using Theorem 3 we find 8 conformations for the correlator in Eq. (38) whose exponents read


and Theorem 6 yields their multiplicities


We check the solution by inspecting


where the RHS follows from and .

For the denominator we get




We again check the consistence


By comparing with [26](ch. 7) we also find a perfect agreement. From this point we proceed as usual. Since we integrate over the ‘internal’ variables (here labelled 3 and 4), some conformations represent the same physical process. We can easily find them from (40) as those invariant wrt the swap . Those are the following pairs of ten-tuples: , and . Considering the double multiplicities we reproduce and finalize the derivation of the second order of the theory [26](ch. 7). ∎

Of course, these are standard results obtainable through Wick’s theorem. Yet, the ease with which we reproduced them is notable as we don’t need to count the contractions or draw the Feynman diagrams. An implementation of Theorems 3 and 6 in a suitable programming language such as Mathematica is straightforward. For example, one can convert the proof of Theorem 3 into a program that systematically finds all solutions. Essentially, if it takes one time step to find the first solution then all solutions can be find in time steps where increases polynomially. There is really no Diophantine calculation involved since due to the simple form of (17) we get all the solutions by inspection and only need to list them (save them to a memory).

The whole process would arguably be much more impressive for more than two interacting fields and for high-loop corrections. In that case the number of variables increases and Theorem 3 needs to be generalized. This can be done in a straightforward way but the presence of more variables is ultimately responsible for the factorial growth of the perturbative contributions. We obviously cannot get rid of the factorial growth but the easy-to-see character of the Diophantine linear system remains and so all the solutions can again be listed in time steps ( grows factorially with a perturbative order, though). In this way, we can again get the perturbative contributions for any number of interacting fields, to a essentially an arbitrarily high perturbative order (keeping in mind the exponential growth discussed in this paragraph) and for all with no effort whatsoever.

Furthermore, the Feynman diagrams are encoded in the conformation exponents. This is not surprising as represents a free propagator for being the internal degrees of freedom. So we easily recover them from our approach if they are still needed. Let’s illustrate it on our next example.

Example 5 ( for