format: text

text:
# Quaternionic Madelung Transformation and Non-Abelian Fluid Dynamics

## I Introduction

Since the $`1930`$’s it has been known that it is possible to construct quantum mechanics over fields other than the complex numbers. In particular, a detailed and consistent presentation of quaternionic quantum mechanics has been developed and found to have a number of fascinating theoretical properties Adler (1995). In the present work, we develop a relationship between quaternionic quantum theory and the simulation of non-relativistic, non-abelian fluid dynamics. Many computer simulation techniques developed for complex quantum mechanics generalize straightforwardly to the quaternionic case, so the observation made in this paper makes possible new simulation strategies for non-abelian fluid dynamics.

The Madelung transformation transforms the complex Schroedinger equation into the continuity and Euler equations describing an inviscid fluid. The magnitude squared of the wavefunction is the mass density, and the gradient of the phase of the wavefunction is the hydrodynamic velocity. The mass continuity equation follows directly from the conservation of probability density, and the momentum equation of motion follows with the identification of a density dependent pressure.

Recently, there has been interest in developing non-abelian fluid dynamical models to describe systems such as quark-gluon plasmas, where one may wish the color degrees of freedom to be retained in the fluid description Jackiw et al. (2000); Bistrovic et al. (2002); Jackiw (2002). The color degrees of freedom are associated with the gauge symmetry group of the field; for the example of a quark-gluon plasma this would be the group $`SU(3)`$ of quantum chromodynamics (QCD).

A natural question to ask in the context of quaternionic quantum mechanics is: What is the generalization of the Madelung transformation for the quaternionic Schroedinger equation? As the quaternionic elements are anticommuting, and in particular have the Lie algebra of $`SU(2)`$, we might plausibly expect the fluid equations arising from such a transformation to be those of a non-abelian fluid.

The remainder of this paper is organized as follows. We first identify the quaternionic generalization of the Madelung transformation, and identify the hydrodynamic variables. In order to find evolution equations for these quantities we cast the quaternionic Schroedinger equation in Hamiltonian form, defining a canonical Poisson bracket and Hamiltonian for the quaternionic wavefunction. We then use our identification of the quaternionic Madelung variables to define a non-canonical transformation of our Hamiltonian field theory, obtaining the bracket structure and evolution equations. We compare our approach with that of Jackiw et al. Jackiw et al. (2000); Bistrovic et al. (2002); Jackiw (2002), and briefly present several numerical algorithms by which our equations may be straightforwardly simulated. We close the paper with some conclusions and directions for future work.

## II Madelung variables

The identification of the gradient of the phase as the relevant hydrodynamic variable in the complex case is a strong indication that local gauge symmetry plays an important role in the transformation, and will therefore be useful as a guide in the more difficult quaternionic case. We first identify the gradient of the phase of the wavefunction, which is identified with the fluid velocity in the complex Madelung transformation, as a pure $`U(1)`$ gauge. The condition that the fluid velocity be a pure gauge translates to the constraint of vanishing vorticity for the Madelung fluid. We then use this observation to obtain the form of a pure gauge for a quaternionic $`SU(2)`$ Yang-Mills gauge field. We find three vector fields and show that they obey constraints which are the generalization of vanishing vorticity to the non-abelian case.

To perform the program outlined above we first construct the gauge transformation of a gauge field $`A_i`$ with arbitrary gauge symmetry group $`𝒢`$. Our gauge field is introduced via the gauge covariant derivative $`D_\alpha =_\alpha A_\alpha `$. All physical quantities are invariant under the transformations

$$\psi \psi ^{}=g\psi $$
(1)
and

$$A_\alpha A_\alpha ^{},$$
(2)
where $`g`$ is an element of the gauge symmetry group and the transformation properties of $`A_i`$ are obtained by requiring

$$\left(_\alpha A_\alpha ^{}\right)g\psi =g\left(_\alpha A_\alpha \right)\psi ,$$
(3)
which gives in turn

$$A_\alpha ^{}=\left(_\alpha g\right)g^1+gA_\alpha g^1.$$
(4)
For complex quantum mechanics coupled to the electromagnetic field the gauge group $`𝒢`$ is $`U(1)`$, the element $`g`$ is any unimodular complex number, and the quantity $`\left(_\mu g\right)g^1`$ is simply the gradient of the change in phase of the wavefunction. The above transformation must produce no change in electromagnetic field strength $`F_{\alpha \beta }=_\alpha A_\beta _\beta A_\alpha `$. It is clear that they do not, as the change in $`A_\mu `$ takes the form of a similarity transformation and the addition of the pure gauge $`\left(_\alpha g\right)g^1`$. The condition that the gauge fields remain unchanged is the elementary condition $`\times \theta =0`$, i.e., the constraint that the Madelung fluid has vanishing vorticity.

The above discussion motivates a straightforward generalization of the Madelung transformation. The Madelung velocity fields have the form of pure $`U(1)`$ gauges, and the condition that pure gauges do not contribute to the gauge field strengths translates into the condition of vanishing vorticity for the Madelung velocity field. Turning to the quaternionic case, the gauge group $`𝒢`$ is $`SU(2)`$, the element $`g`$ is any unimodular quaternion, and we identify the pure gauges $`\left(_\alpha g\right)g^1`$ with the relevant hydrodynamic variables for the quaternionic Madelung transformation.

We write the group element $`g`$ as

$$g=e^{i\mu }e^{j\theta }e^{k\nu }$$
(5)
Where $`i`$, $`j`$ and $`k`$ are the quaternionic elements satisfying $`i^2=j^2=k^2=ijk=1`$, and $`\mu `$, $`\nu `$ and $`\theta `$ are real angles. Substitution of this form for $`g`$ into Eq. (5) gives

$$(g)g^1=i\mu +e^{i\mu }je^{i\mu }\theta +e^{i\mu }e^{j\theta }ke^{j\theta }e^{i\mu }\nu ,$$
(6)
We write

$$(g)g^1=i𝐮+j𝐯+k𝐰,$$
(7)
where

$$\begin{array}{cc}\hfill 𝐮& =\mu +\mathrm{cos}2\theta \nu \hfill \\ \hfill 𝐯& =\mathrm{cos}2\mu \theta +\mathrm{sin}2\theta \mathrm{sin}2\mu \nu \hfill \\ \hfill 𝐰& =\mathrm{sin}2\mu \theta \mathrm{sin}2\theta \mathrm{cos}2\mu \nu .\hfill \end{array}$$
(8)
The definition of the density is $`\rho =\psi \overline{\psi }`$, where $`\psi `$ is the quaternionic wavefunction $`\psi =\sqrt{\rho }g`$ and its quaternionic conjugate $`\overline{\psi }`$ is $`\psi =\sqrt{\rho }g^1=\sqrt{\rho }e^{k\nu }e^{j\theta }e^{i\mu }`$. The quaternionic generalizations of the Madelung variables are then $`\rho `$, $`𝐮`$, $`𝐯`$ and $`𝐰`$.

These quaternionic Madelung variables are the pure gauges of an $`SU(2)`$ Yang-Mills theory. We can see this by showing that the Yang-Mills field strengths are zero. The non-abelian generalization of the Maxwell field strength is

$$F_{\alpha \beta }=_\alpha A_\beta _\beta A_\alpha [A_\alpha ,A_\beta ]$$
(9)
where $`A_\alpha `$ are quaternionic gauge fields. Sufficient and necessary conditions for these field strengths to be zero are

$$\begin{array}{cc}\hfill \times 𝐮& =2𝐯\times 𝐰\hfill \\ \hfill \times 𝐯& =2𝐰\times 𝐮\hfill \\ \hfill \times 𝐰& =2𝐮\times 𝐯\hfill \end{array}$$
(10)
It may be readily verified that these equations are satisfied by $`𝐮`$, $`𝐯`$ and $`𝐰`$, and hence are the quaternionic generalization of vanishing vorticity in the complex Madelung transformation. We now obtain the equations of motion for $`𝐮`$, $`𝐯`$ and $`𝐰`$.

## III Equations of Motion

In order to obtain the equations of motion for $`𝐮,𝐯`$ and $`𝐰`$ we first obtain a canonical Hamiltonian form for the quaternionic Schroedinger equation. To do this we need to define a Hamiltonian and a Poisson bracket. Defining the Hamiltonian is straightforward, as it follows directly from the definition for the complex case,

$$\begin{array}{cc}\hfill H& =\frac{1}{2}\overline{\psi }\widehat{H}\psi d^3z=\frac{1}{2}\left(\frac{1}{2}\overline{\psi }^2\psi +V\overline{\psi }\psi \right)d^3z\hfill \\ & =\frac{1}{2}\left(\frac{|\rho |^2}{8\rho }+\frac{\rho }{2}(𝐮^2+𝐯^2+𝐰^2)+V\rho \right)d^3x,\hfill \end{array}$$
(11)
where we use units such that $`\mathrm{}=m=1`$ throughout. By the argument presented in Appendix A we define the bracket

$$\begin{array}{cc}\hfill \{Q(𝐱),P(𝐲)\}& =\left(\frac{\delta Q(𝐱)}{\delta a(𝐳)}\frac{\delta P(𝐲)}{\delta b(𝐳)}\frac{\delta Q(𝐱)}{\delta b(𝐳)}\frac{\delta P(𝐲)}{\delta a(𝐳)}\right)d^3z\hfill \\ & +\left(\frac{\delta Q(𝐱)}{\delta d(𝐳)}\frac{\delta P(𝐲)}{\delta c(𝐳)}\frac{\delta Q(𝐱)}{\delta c(𝐳)}\frac{\delta P(𝐲)}{\delta d(𝐳)}\right)d^3z,\hfill \end{array}$$
(12)
where $`a`$, $`b`$ , $`c`$ and $`d`$ are the quaternionic components of the wavefunction $`\psi =a+ib+jc+kd`$. The above Poisson bracket and Hamiltonian reproduce the quaternionic Schroedinger equation. The definitions Eq. (8) are then a non-canonical transformation of this Hamiltonian field theory. The bracket structure of the quaternionic Madelung equations follows by lengthy but straightforward algebra. We spare the reader the details and present the noncanonical brackets amongst $`\rho `$, $`𝐮`$, $`𝐯`$ and $`𝐰`$ in Appendix B. Using the above Hamiltonian, the equation of motion for $`\rho `$ is then

$$\frac{\rho }{t}+\left(\rho 𝐮\right)=0.$$
(13)
This is the continuity equation, and if $`\rho `$ is identified as the singlet density, this equation identifies $`𝐮`$ as the singlet velocity. The equation of motion for $`𝐮`$ is

$$\begin{array}{cc}\hfill \frac{𝐮}{t}& =V^{}\frac{1}{2}(𝐮^2+𝐯^2+𝐰^2)\frac{𝐯}{\rho }(\rho 𝐯)\frac{𝐰}{\rho }(\rho 𝐰),\hfill \end{array}$$
(14)
where $`V^{}`$ is the density-dependent potential (also called the Bohmian potential) which also appears in the complex Madelung transformation,

$$V^{}=V\frac{1}{4\sqrt{\rho }}\left(\frac{\rho }{\sqrt{\rho }}\right).$$
(15)
Note that the curl of Eq. (14) is generally nonvanishing. Since $`𝐮`$ is identified as the singlet density, this means that the quaternionic Madelung fluid has non-zero vorticity. This can be clearly seen by inspecting the $`\{𝐮(𝐱),𝐮(𝐲)\}`$ bracket, which has the correct form for rotational flow. The equations for $`𝐯`$ and $`𝐰`$ are

$$\begin{array}{cc}\hfill \frac{𝐯}{t}& =\frac{1}{2}\left[\frac{(\rho 𝐰)}{\rho }\right]2𝐰\left[V^{}(𝐱)+\frac{1}{2}(𝐮^2+𝐯^2+𝐰^2)\right]+\frac{𝐮}{\rho }(\rho 𝐰)\hfill \end{array}$$
(16)
and

$$\begin{array}{cc}\hfill \frac{𝐰}{t}& =\frac{1}{2}\left[\frac{(\rho 𝐯)}{\rho }\right]+2𝐯\left[V^{}(𝐱)+\frac{1}{2}(𝐮^2+𝐯^2+𝐰^2)\right]+\frac{𝐮}{\rho }(\rho 𝐯).\hfill \end{array}$$
(17)
One immediately notices that the equations for $`𝐯`$ and $`𝐰`$ are invariant under the duality transformation

$$\left(\begin{array}{c}𝐯^{}\\ 𝐰^{}\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}\theta & \mathrm{sin}\theta \\ \mathrm{sin}\theta & \mathrm{cos}\theta \end{array}\right)\left(\begin{array}{c}𝐯\\ 𝐰\end{array}\right)$$
(18)
This duality symmetry is a consequence of the global gauge invariance of the quaternionic Schroedinger equation. This global invariance is discussed in detail in Appendix C.

## IV Non-abelian fluid dynamics

Recently there has been interest in constructing non-abelian fluid dynamics and fluid magnetohydrodynamics in order to simulate the long-wavelength degrees of freedom in a quark-gluon plasma. In particular, Jackiw et al. have described several models for the relativistic dynamics of such a system Jackiw et al. (2000); Jackiw (2002). Those authors have also presented a non-relativistic multicomponent wavefunction very similar to the quaternionic wavefunction considered here Bistrovic et al. (2002). In fact, as we shall show, their model differs from the quaternionic Madelung transformation we have investigated only by a global gauge of the type described in detail in Appendix C. In this section we compare their multicomponent wave equation with our approach and give the global gauge transformation relating their work to our own.

Bistrovic et al. consider a two component complex wave equation,

$$i\frac{d\psi }{dt}=\frac{1}{2}^2\psi .$$
(19)
Those authors consider the group $`SU(2)`$ and the fundamental representation: $`T^a=\sigma ^a/(2i)`$, $`\{T^a,T^b\}=\delta ^{ab}/2`$, where $`\sigma ^a,a=1,2,3`$ are the Pauli matrices. The $`T^a`$’s here are simply a two dimensional complex representation of the quaternionic elements

$$T^1=\frac{i}{2},T^2=\frac{j}{2},T^3=\frac{k}{2}.$$
(20)
Those authors define the conserved singlet and color current;

$$j^\mu =(c\rho ,𝐣),\rho =\psi ^{}\psi ,𝐣=\mathrm{Im}\psi ^{}\psi .$$
(21)
$$J_a^\mu =(c\rho _a,𝐉_a),\rho _a=i\psi ^{}T^a\psi ,𝐉_a=\mathrm{Re}\psi ^{}T^a\psi .$$
(22)
and the wavefunction,

$$\psi =\sqrt{\rho }g𝐀$$
(23)
where $`\rho `$ is the scalar $`\psi ^{}\psi `$, where $``$ denotes hermitian conjugation, $`g`$ is a group element, and $`𝐀`$ is an arbitrary constant vector taken to be

$$𝐀=\left(\begin{array}{c}1\\ 0\end{array}\right)$$
(24)
in which case $`i𝐀^{}T^a𝐀=\delta ^{a3}/2`$. The singlet density is $`\rho `$, while the singlet current $`𝐣`$ is

$$𝐣=\rho 𝐩,𝐩i𝐀^{}g^1g𝐀.$$
(25)
The key relation connecting our approach to Jackiw’s is:

$`g^1g`$ $``$ $`2𝐩^aT^a`$ (26)
$`𝐩`$ $`=`$ $`𝐩^at^a,`$ (27)
where $`t^a/2=iu^{}T^au=\delta ^{a3}/2`$, and hence

$$𝐩=𝐩^at^a=𝐩^3.$$
(28)
It follows that the singlet velocity in Bistrovic et al. (2002) is equated to the negative of third quaternionic component of the wavefunction in our approach. His multicomponent wavefunction is therefore equivalent to a quaternionic Schroedinger equation

$$k\frac{d\psi }{dt}=\frac{1}{2}^2\psi .$$
(29)
This quaternionic Schroedinger equation is related to the one considered in this paper by a one-parameter family of global gauge transformations, described in detail in Appendix C.

The color and singlet currents identified in Bistrovic et al. (2002) may be obtained from the fields $`𝐮`$, $`𝐯`$, and $`𝐰`$ defined for our quaternionic Madelung equation. The quaternionic Madelung transformation we have described is a particularly simple implementation of this type of non-abelian fluid model, but its usefulness is limited as the approach cannot be extended to groups other than $`SU(2)`$. The quaternionic approach has some advantages in terms of ease of numerical simulation, which we describe below.

## V Simulation

The Schroedinger equation may be simulated by a variety of methods. We briefly describe three methods, all of which may be trivially generalized to the simulation of quaternionic quantum mechanics. The first method is operator splitting, or its higher order variant, the symplectic integrator method. This method is perfectly unitary, can be generalized to three dimensions, and can be generalized to high orders of accuracy Chin and Chen (2001); Feit et al. (1982a, b); Takahashi and Ikeda (1993); Gray and Manolopoulos (1996); Takahashi and Ikeda (1997). The second-order spatial derivatives required in the operator splitting method may be avoided by using a quantum lattice-Boltzmann method (QLBE) Succi and Benzi (1993); Succi (1996, 2002). The QLBE is unconditionally stable, perfectly unitary, and offers considerable advantages for parallel computing applications.

Closely related to the quantum lattice Boltzmann approach is the quantum lattice-gas model. These algorithms were originally developed for simulation of the Dirac equation in one spatial dimension Meyer (1996a, b, 1997a, 1997b). The non-relativistic quantum lattice gas was subsequently developed for the simulation of many-body quantum mechanics in arbitrary numbers of dimensions Boghosian and Taylor (1998). The quantum lattice-gas algorithm for a single nonrelativistic particle provides a tractable technique for the simulation of the quaternionic Schroedinger equation.

It should be noted that the quantum lattice gas algorithm (QLGA) of Boghosian and Taylor (1998) was developed in order to find a quantum algorithm for Schroedinger’s equation. The QLGA provides a method of simulating quantum mechanics on a quantum computer with an exponential enhancement in the performance. Minor modifications to this QLGA should also provide a method of simulating the quaternionic Madelung equations derived above on a (complex) quantum computer with a similar exponential speedup.

## VI Conclusions

We have obtained the quaternionic generalization of the Madelung transformation. We have identified the generalizations of the hydrodynamic fields and obtained their evolution equations by a non-canonical transformation of the Hamiltonian field theory for the quaternionic Schroedinger equation.

The equations we have obtained also possess a limit which may be of interest. The choice of arbitrarily large potential may always make the second terms in Eqs. (16) and (17) dominant. In this case a solution which is constant in space and oscillatory in time will be a good first-order approximation to the behavior of $`𝐯`$ and $`𝐰`$. This suggests that there may exist an approximate treatment of these equations in which $`𝐯`$ and $`𝐰`$ are treated as “fast” variables which may be averaged away. This would be analogous to the gyro-averaged description of MHD, in which the fast gyrofrequency motions are averaged away and only the slow motion of the gyro centers is retained in the reduced dynamics Littlejohn (1983, 1979); Boghosian (1987). The Hamiltonian approach taken in this paper provides a good starting point for the derivation of such an averaged treatment.

The Madelung transformation of complex quantum mechanics seems intially computationally promising, as it shows us how to simulate the nonlinear Euler equations by simulating the linear Schroedinger equation. However, irrotational inviscid fluid flow is not a problem of sufficient interest to merit serious computational investigation. The quaternionic Madelung transformation gives rise to equations of considerably more interest, as they describe the (rotational) flow of a non-abelian fluid. However, the constraints in Eq. (10) are the restriction on the non-abelian fluid, analogous to vanishing vorticity in the abelian case. This suggests that there exists a generalization of the equations obtained here to non-abelian flow fields which do not satisfy the constraints of Eq. (10). It is possible that these general non-abelian flow equations will contain as much additional physics as does rotational Euler flow compared with irrotational Euler flow.

A natural question arises as to whether one may generalize this approach to groups other than $`SU(2)`$ and fields other than the quaternions. We may answer this immediately in the negative: There are no division algebras other than the reals, the complex numbers, the quaternions and the the octonions Kantor and Solodovnikov (1973). The octonions form a division algebra, but are not associative, and this precludes the formulation of an octonionic quantum mechanics Adler (1995).

Finally, we have noted that the computational tractability of the complex Schroedinger equation is shared by its quaternionic generalization, and so we anticipate that extensive simulation of the quaternionic Schroedinger equation by any one of the methods described above will yield results of interest to the particle physics community.

## Acknowledgements

BMB was supported in part by the U.S. Air Force Office of Scientific Research under grant number F49620-01-1-0385. PJL was supported by the DARPA QuIST program under AFOSR grant number F49620-01-1-0566.

## Appendix A Quaternionic Poisson Bracket

If we write the Schroedinger equation as a pair of real equations for the real and imaginary parts of the wavefunction $`\psi =a+ib`$ we obtain:

$$\frac{d}{dt}\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{cc}0& \widehat{H}\\ \widehat{H}& 0\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right)$$
(30)
where $`\widehat{H}=\frac{1}{2}^2+V`$. Note that the above construction does not enable us to construct quantum mechanics over the reals. The symplectic matrix in the above equation does not have real eigenvalues and so there are no energy eigenstates in real quantum mechanics Adler (1995). However, the form of this equation makes it immediately obvious that the Schroedinger equation may be described by a canonical Hamiltonian field theory in which the real and imaginary parts of the equation are canonically conjugate. The definition of the bracket follows automatically,

$$\begin{array}{cc}\hfill \{Q(𝐱),P(𝐲)\}& =\left(\frac{\delta Q(𝐱)}{\delta a(𝐳)}\frac{\delta P(𝐱)}{\delta b(𝐳)}\frac{\delta Q(𝐱)}{\delta b(𝐳)}\frac{\delta P(𝐱)}{\delta a(𝐳)}\right)d^3z\hfill \end{array}$$
(31)
The quaternionics have a four dimensional real representation, which enables us to write the quaternionic Schroedinger equation as four real equations for the four components of the wavefunction $`\psi =a+ib+jc+kd`$,

$$\frac{d}{dt}\left(\begin{array}{c}a\\ b\\ c\\ d\end{array}\right)=\left(\begin{array}{cccc}0& \widehat{H}& 0& 0\\ \widehat{H}& 0& 0& 0\\ 0& 0& 0& \widehat{H}\\ 0& 0& \widehat{H}& 0\end{array}\right)\left(\begin{array}{c}a\\ b\\ c\\ d\end{array}\right)$$
(32)
Again, we see that these equations have a symplectic structure where $`a`$ and $`b`$ are canonically conjugate, and so are $`d`$ and $`c`$. The definition of the bracket, Eq. (12), again follows automatically.

## Appendix B Poisson Bracket structure

Our quaternionic Madelung variables have the following bracket structure:

$$\begin{array}{cc}\hfill \{\rho (𝐱),w_i(𝐲)\}& =4\delta (𝐱𝐲)v_i(𝐱)\hfill \\ \hfill \{\rho (𝐱),v_i(𝐲)\}& =4\delta (𝐱𝐲)w_i(𝐱)\hfill \\ \hfill \{\rho (𝐱),u_i(𝐲)\}& =2\delta _i^{}(𝐱𝐲)\hfill \end{array}$$
(33)
$$\begin{array}{cc}\hfill \{u_i(𝐱),u_j(𝐲)\}& =\frac{2\delta (𝐱𝐲)}{\rho (𝐱)}\left(_iu_j(𝐱)_ju_i(𝐱)\right)\hfill \\ \hfill \{u_i(𝐱),v_j(𝐲)\}& =\frac{2\delta _j^{}(𝐱𝐲)v_i(𝐱)}{\rho (𝐱)}+\frac{4\delta (𝐱𝐲)}{\rho (𝐱)}\left[\frac{1}{2}\left(_iv_j(𝐱)_jv_i(𝐱)\right)u_i(𝐱)w_j(𝐱)\right]\hfill \\ \hfill \{u_i(𝐱),w_j(𝐲)\}& =\frac{2\delta _j^{}(𝐱𝐲)w_i(𝐱)}{\rho (𝐱)}\frac{4\delta (𝐱𝐲)}{\rho (𝐱)}\left[\frac{1}{2}(_jw_i(𝐱)_iw_j(𝐱))v_j(𝐱)u_i(𝐱)\right]\hfill \\ \hfill \{v_i(𝐱),u_j(𝐲)\}& =\frac{2\delta _i^{}(𝐲𝐱)v_j(𝐲)}{\rho (𝐲)}\frac{4\delta (𝐲𝐱)}{\rho (𝐲)}\left[\frac{1}{2}\left(_jv_i(𝐲)_iv_j(𝐲)\right)u_j(𝐲)w_i(𝐲)\right]\hfill \\ \hfill \{v_i(𝐱),v_j(𝐲)\}& =\frac{2\delta _i^{}(𝐲𝐱)u_j(𝐲)}{\rho (𝐲)}+\frac{2\delta _j^{}(𝐱𝐲)u_i(𝐱)}{\rho (𝐱)}\hfill \\ & \frac{4\delta (𝐲𝐱)}{\rho (𝐱)}\left[\frac{1}{2}\left(_iu_j(𝐱)_ju_i(𝐱)\right)+v_i(𝐱)w_j(𝐱)v_j(𝐱)w_i(𝐱)\right]\hfill \\ \hfill \{v_i(𝐱),w_j(𝐲)\}& =\frac{\delta _{ij}^{\prime \prime }(𝐱𝐲)}{\rho (𝐱)}\frac{\delta _j^{}(𝐱𝐲)_i\rho (𝐱)}{\rho (𝐱)^2}4\delta (𝐲𝐱)\frac{u_j(𝐱)u_i(𝐱)}{\rho (𝐱)}\hfill \\ \hfill \{w_i(𝐱),u_j(𝐲)\}& =\frac{2\delta _i^{}(𝐲𝐱)w_j(𝐲)}{\rho (𝐲)}+\frac{4\delta (𝐱𝐲)}{\rho (𝐱)}\left[\frac{1}{2}(_iw_j(𝐱)_jw_i(𝐱))v_i(𝐱)u_j(𝐱)\right]\hfill \\ \hfill \{w_i(𝐱),v_j(𝐲)\}& =\frac{\delta _{ij}^{\prime \prime }(𝐲𝐱)}{\rho (𝐲)}+\frac{\delta _i^{}(𝐲𝐱)_j\rho (𝐲)}{\rho (𝐲)^2}+4\delta (𝐱𝐲)\frac{u_i(𝐱)u_j(𝐱)}{\rho (𝐱)}\hfill \\ \hfill \{w_i(𝐱),w_j(𝐲)\}& =2\frac{\delta _i^{}(𝐲𝐱)}{\rho (𝐲)}u_j(𝐲)+2\frac{\delta _j^{}(𝐱𝐲)}{\rho (𝐱)}u_i(𝐱)\hfill \\ & +4\delta (𝐲𝐱)\frac{1}{\rho (𝐲)}\left[\frac{1}{2}(_iu_j(𝐲)_ju_i(𝐲))+v_i(𝐲)w_j(𝐲)v_j(𝐲)w_i(𝐲)\right]\hfill \end{array}$$
(34)
## Appendix C Global gauge symmetry

We obtained our quaternionic Madelung variables by considering the local gauge invariance properties of our Schroedinger equation coupled to a quaternionic gauge field. The observables of quaternionic quantum mechanics posess global $`SU(2)`$ gauge invariance, and this invariance will lead to a set of symmetry properties for the quaternionic Madelung equations. The Schroedinger equation is

$$i\frac{\psi }{t}=\widehat{H}\psi $$
(35)
Under the global gauge transformation $`\psi \psi ^{}=\varphi \psi `$ this becomes:

$$\varphi ^1i\varphi \frac{\psi }{t}=\widehat{H}\psi $$
(36)
If $`\varphi `$ is a unimodular quaternion then the transformation $`\varphi ^1Q\varphi `$, where $`Q`$ is any quaternion, is an automorphism of the quaternions, and $`\varphi ^1Q\varphi `$ is also a unimodular quaternion. The equivalent global gauge transformation applied to the complex Schroedinger equation also induces an automorphism of the complex numbers, but in that case it is the trivial automorphism $`\theta :e^{i\theta }Ze^{i\theta }=Z`$. The choice of $`i`$ in Eq. (35), which we inherited from the complex Schroedinger equation, constitutes a global gauge choice. We could equally well write the quaternionic Schroedinger equation replacing $`i`$ by $`j`$, or $`k`$, or indeed any pure imaginary unimodular quaternion. We now consider how such global gauge choices affect our quaternionic Madelung equations.

Our Madelung variables are defined by:

$$(g)g^1=i𝐮+j𝐯+k𝐰,$$
(37)
applying a global gauge transformation to this quantity gives:

$$\varphi ^1(g)g^1\varphi =\varphi ^1i\varphi 𝐮+\varphi ^1j\varphi 𝐯+\varphi ^1k\varphi 𝐰,$$
(38)
Taking the quaternion automorphism defined by:

$$f_p=\overline{\varphi }e_p\varphi .$$
(39)
Where $`e_p=i,j,k,p=1,2,3`$ are our original quaternionic elements and $`f_q`$ are our new quaternionic elements. We express the new quaternionic elements in terms of the old:

$$f_q=a_{qp}e_p$$
(40)
We can determine the properties of the coefficients $`a_{qp}`$ from the algebra obeyed by both sets of quaternionic elements:

$$\begin{array}{cc}\hfill e_qe_r& =\delta _{qr}+ϵ_{qrs}e_s\hfill \end{array}$$
(41)
Substituting:

$$\begin{array}{cc}\hfill a_{mp}a_{nq}\left(\delta _{pq}+ϵ_{pqs}e_s\right)& =\delta _{mn}+ϵ_{mno}a_{or}e_r\hfill \end{array}$$
(42)
Where $`ϵ_{mno}`$ is the completely antisymmetric Levi-Civita tensor. Separating the coefficients of the quaternionic elements $`e_s`$ gives:

$$\begin{array}{cc}\hfill a_{rp}a_{sq}\delta _{pq}& =\delta _{rs}\hfill \\ \hfill (a_{rp}a_{sq}ϵ_{pqr}& ϵ_{rsl}a_{lr})e_m=0\hfill \end{array}$$
(43)
These equations show that the matrix $`a_{pq}`$ is orthogonal with unit determinant. The global gauge transformation induces a linear transformation (actually a proper rotation) of our Madelung variables,

$$\begin{array}{cc}\hfill 𝐮^{}& =a_{11}𝐮+a_{21}𝐯+a_{31}𝐰\hfill \\ \hfill 𝐯^{}& =a_{12}𝐮+a_{22}𝐯+a_{32}𝐰\hfill \\ \hfill 𝐰^{}& =a_{13}𝐮+a_{23}𝐯+a_{33}𝐰.\hfill \end{array}$$
(44)
The Madelung transformation therefore applies to a three-parameter family of quaternionic Schroedinger equations, and the results of simulating any member of this family may be related to the results of simulating any other member by a global gauge transformation. As a specific example of this we construct the transformation relating the quaternionic Schroedinger equations:

$$i\frac{d\psi }{dt}=\frac{1}{2}^2\psi $$
(45)
and

$$k\frac{d\psi }{dt}=\frac{1}{2}^2\psi .$$
(46)
The two equations are related by the global gauge transformation $`\psi ^{}\varphi \psi `$ where $`\varphi =e^{i\theta }\left(1j\right)/\sqrt{2}`$, and the matrix $`a_{pq}`$ is

$$a_{pq}=\left(\begin{array}{ccc}0& 0& 1\\ 0& 1& 0\\ 1& 0& 0\end{array}\right)$$
(47)
This is the global gauge transformation relating the quaternionic Schroedinger equation considered here to the multicomponent non-abelian Schroedinger equation considered by Bistrovic et al. Bistrovic et al. (2002).

The gauge transformations relating the two quaternionic Schroedinger equations above form a one parameter family. There is always a one-parameter subgroup of the full gauge group which leaves any quaternionic Schroediger equation identically invariant. For the Schroedinger equation, Eq. (45), the elements of this group may be written $`g=e^{i\theta /2}`$. This subgroup induces the transformations

$$a_{pq}=\left(\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}\theta & \mathrm{sin}\theta \\ 0& \mathrm{sin}\theta & \mathrm{cos}\theta \end{array}\right)$$
(48)
This is the set of global gauge transformations which give rise to the duality symmetry of $`𝐯`$ and $`𝐰`$ noted in Eq. (18).


text:
# 1 Introduction

## 1 Introduction

The strict positive realness (SPR) of transfer functions is an important performance specification, and plays a critical role in various fields such as absolute stability/hyperstability theory , passivity analysis , quadratic optimal control  and adaptive system theory . In recent years, stimulated by the parametrization method in robust stability analysis , the study of robust strictly positive real systems has received much attention, and great progress has beem made  - - - -. However, most results belong to the category of robust SPR analysis. Valuable results in robust SPR synthesis are rare. The following fundamental problem is still open  -  -:

Suppose $`a(s)`$ and $`b(s)`$ are two $`n`$-$`th`$ order Hurwitz polynomials, does there exist, and how to find a (fixed) polynomial $`c(s)`$ such that $`c(s)/a(s)`$ and $`c(s)/b(s)`$ are both SPR?

By the definition of SPR, it is easy to know that the Hurwitz stability of the convex combination of $`a(s)`$ and $`b(s)`$ is necessary for the existence of polynomial $`c(s)`$ such that $`c(s)/a(s)`$ and $`c(s)/b(s)`$ are both SPR. In , it was proved that, if $`a(s)`$ and $`b(s)`$ have the same even (or odd) parts, such a polynomial $`c(s)`$ always exists; In , it was proved that, if $`n4`$ and $`a(s),b(s)K`$ ($`K`$ is a stable interval polynomial set), such a polynomial $`c(s)`$ always exists; Recent results show that  -, if $`n6`$ and $`a(s)`$ and $`b(s)`$ are the two endpoints of the convex combination of stable polynomials, such a polynomial $`c(s)`$ always exists. Some sufficient condition for robust SPR synthesis are presented in , especially, the design method in  is numerically efficient for high-order polynomial segments or interval polynomials, and the derived conditions are necessary and sufficient for low-order polynomial segments or interval polynomials.

This paper presents a constructive proof of the following statement: for any two $`n`$-$`th`$ order polynomials $`a(s)`$ and $`b(s),`$ the Hurwitz stability of their convex combination is necessary and sufficient for the existence of a polynomial $`c(s)`$ such that $`c(s)/a(s)`$ and $`c(s)/b(s)`$ are both SPR. This also shows that the conditions given in  are also necessary and sufficient and the above open problem admits a positive answer. Some previously obtained SPR synthesis results for low-order polynomial segments  - become special cases of our main result in this paper.

## 2 Main Results

In this paper, $`P^n`$ stands for the set of $`n`$-$`th`$ order polynomials of $`s`$ with real coefficients, $`R`$ stands for the field of real numbers, $`(p)`$ stands for the order of polynomial $`p()`$, and $`H^nP^n`$ stands for the set of $`n`$-$`th`$ order Hurwitz stable polynomials with real coefficients.

In the following definition, $`p()P^m,q()P^n,f(s)=p(s)/q(s)`$ is a rational function.

Definition 1  $`f(s)`$ is said to be strictly positive real(SPR), if

(i) $`(p)=(q);`$

(ii) $`f(s)`$ is analytic in $`\text{Re}[s]0,`$ (namely, $`q()H^n`$ );

(iii) $`\text{Re}[f(j\omega )]>0,\omega R.`$

If $`f(s)=p(s)/q(s)`$ is proper, it is easy to get the following property:

Property 1  If $`f(s)=p(s)/q(s)`$ is a proper rational function, $`q(s)H^n,`$ and $`\omega R,\text{Re}[f(j\omega )]>0,`$ then $`p(s)H^nH^{n1}.`$

The following theorem is the main result of this paper:

Theorem 1 Suppose $`a(s)=s^n+a_1s^{n1}+\mathrm{}+a_nH^n,b(s)=s^n+b_1s^{n1}+\mathrm{}+b_nH^n,`$ the necessary and sufficient condition for the existence of a polynomial $`c(s)`$ such that $`c(s)/a(s)`$ and $`c(s)/b(s)`$ are both Strict Positive Real, is

$$\lambda b(s)+(1\lambda )a(s)H^n,\lambda [0,1].$$

The statement is obviously true for the cases when $`n=1`$ or $`n=2.`$ We will prove it for the case when $`n3.`$

Since SPR transfer functions enjoy convexity property, by Property 1, we get the necessary part of the theorem.

To prove sufficiency, we must first introduce some lemmas.

Lemma 1 Suppose $`a(s)=s^n+a_1s^{n1}+\mathrm{}+a_nH^n,`$ then, for every $`k\{1,2,\mathrm{},n2\},`$ the following quadratic curve is an ellipse in the first quadrant (i.e., $`x_i0,i=1,2,\mathrm{},n1`$) of the $`R^{n1}`$ space $`(x_1,x_2,\mathrm{},x_{n1})`$ 222 When $`n=3,`$ the ellipse equation is (see  for details):
$$(a_2x_1a_1x_2a_3)^24(a_1x_1)(a_3x_2)=0.$$
When $`n=4,`$ the two ellipse equations are (see  - for details):
$$\begin{array}{c}\{\begin{array}{c}(a_2x_1+x_3a_1x_2a_3)^24(a_1x_1)(a_3x_2a_2x_3a_4x_1)=0,\hfill \\ a_4x_3=0,\hfill \end{array}\hfill \\ \{\begin{array}{c}(a_3x_2a_2x_3a_4x_1)^24(a_2x_1+x_3a_1x_2a_3)a_4x_3=0,\hfill \\ a_1x_1=0.\hfill \end{array}\hfill \end{array}$$
When $`n=5,`$ the three ellipse equations are (see  for details):
$$\begin{array}{c}\{\begin{array}{c}(a_2x_1+x_3a_1x_2a_3)^24(a_1x_1)(a_5+a_3x_2+a_1x_4a_2x_3a_4x_1)=0,\hfill \\ a_4x_3a_3x_4a_5x_2=0,a_5x_4=0,\hfill \end{array}\hfill \\ \{\begin{array}{c}(a_5+a_3x_2+a_1x_4a_2x_3a_4x_1)^24(a_2x_1+x_3a_1x_2a_3)(a_4x_3a_3x_4a_5x_2)=0,\hfill \\ a_1x_1=0,a_5x_4=0,\hfill \end{array}\hfill \\ \{\begin{array}{c}(a_4x_3a_3x_4a_5x_2)^24(a_5+a_3x_2+a_1x_4a_2x_3a_4x_1)a_5x_4=0,\hfill \\ a_1x_1=0,a_2x_1+x_3a_1x_2a_3=0.\hfill \end{array}\hfill \end{array}$$
When $`n=6,`$ the four ellipse equations are (see  for details):
$$\begin{array}{c}\{\begin{array}{c}(a_2x_1+x_3a_1x_2a_3)^24(a_1x_1)(a_5+a_3x_2+a_1x_4x_5a_2x_3a_4x_1)=0,\hfill \\ a_6x_1+a_4x_3+a_2x_5a_3x_4a_5x_2=0,a_5x_4a_4x_5a_6x_3=0,a_6x_5=0,\hfill \end{array}\hfill \\ \{\begin{array}{c}(a_5+a_3x_2+a_1x_4x_5a_2x_3a_4x_1)^24(a_2x_1+x_3a_1x_2a_3)(a_6x_1+a_4x_3+a_2x_5a_3x_4a_5x_2)=0,\hfill \\ a_1x_1=0,a_5x_4a_4x_5a_6x_3=0,a_6x_5=0,\hfill \end{array}\hfill \\ \{\begin{array}{c}(a_6x_1+a_4x_3+a_2x_5a_3x_4a_5x_2)^24(a_5+a_3x_2+a_1x_4x_5a_2x_3a_4x_1)(a_5x_4a_4x_5a_6x_3)=0\hfill \\ a_1x_1=0,a_2x_1+x_3a_1x_2a_3=0,a_6x_5=0,\hfill \end{array}\hfill \\ \{\begin{array}{c}(a_5x_4a_4x_5a_6x_3)^24(a_6x_1+a_4x_3+a_2x_5a_3x_4a_5x_2)a_6x_5=0,\hfill \\ a_1x_1=0,a_2x_1+x_3a_1x_2a_3=0,a_5+a_3x_2+a_1x_4x_5a_2x_3a_4x_1=0.\hfill \end{array}\hfill \end{array}$$
:

$`\{\begin{array}{c}c_{k+1}^24c_kc_{k+2}=0,\hfill \\ c_l=0,\hfill \\ l\{1,2,\mathrm{},n\},lk,k+1,k+2,\hfill \end{array}`$

and this ellipse is tangent with the line

$`\{\begin{array}{c}c_l=0,\hfill \\ l\{1,2,\mathrm{},n\},lk+1,k+2,\hfill \end{array}`$

and the line

$`\{\begin{array}{c}c_l=0,\hfill \\ l\{1,2,\mathrm{},n\},lk,k+1,\hfill \end{array}`$

respectively, where $`c_l:=\underset{j=0}{\overset{n}{}}(1)^{l+j}a_jx_{2lj1},l=1,2,\mathrm{},n,`$ $`a_0=1,x_0=1,`$ $`a_i=0`$ if $`i<0`$ or $`i>n,`$ and $`x_i=0`$ if $`i<0`$ or $`i>n1.`$

Proof Since $`a(s)`$ is Hurwitz stable, by using mathematical induction, Lemma 1 is proved by a direct calculation.

For notational simplicity, for $`a(s)=s^n+a_1s^{n1}+\mathrm{}+a_nH^n,b(s)=s^n+b_1s^{n1}+\mathrm{}+b_nH^n,`$ $`k\{1,2,\mathrm{},n2\},`$ denote

$$\begin{array}{cc}\mathrm{\Omega }_{ek}^a:=\{(x_1,x_2,\mathrm{},x_{n1})|\hfill & c_{k+1}^24c_kc_{k+2}<0,\hfill \\ & c_l=0,l\{1,2,\mathrm{},n\},lk,k+1,k+2\},\hfill \end{array}$$

and

$$\begin{array}{cc}\mathrm{\Omega }_{ek}^b:=\{(x_1,x_2,\mathrm{},x_{n1})|\hfill & d_{k+1}^24d_kd_{k+2}<0,\hfill \\ & d_l=0,l\{1,2,\mathrm{},n\},lk,k+1,k+2\},\hfill \end{array}$$

where $`c_l:=\underset{j=0}{\overset{n}{}}(1)^{l+j}a_jx_{2lj1},d_l:=\underset{j=0}{\overset{n}{}}(1)^{l+j}b_jx_{2lj1},l=1,2,\mathrm{},n,`$ $`a_0=1,b_0=1,x_0=1,`$ $`a_i=0`$ and $`b_i=0`$ if $`i<0`$ or $`i>n,`$ and $`x_i=0`$ if $`i<0`$ or $`i>n1.`$

In what follows, $`(A,B)`$ stands for the set of points in the line segment connecting the point $`A`$ and the point $`B`$ in the $`R^{n1}`$ space $`(x_1,x_2,\mathrm{},x_{n1}),`$ not including the endpoints $`A`$ and $`B`$. Denote

$$\begin{array}{cc}\mathrm{\Omega }^a:=\{(x_1,x_2,\mathrm{},x_{n1})|\hfill & (x_1,x_2,\mathrm{},x_{n1})_{i=1,i<jn2}^{n3}(A_i,A_j),\hfill \\ & A_i\mathrm{\Omega }_{ei}^a,i\{1,2,\mathrm{},n2\}\}\hfill \end{array}$$

and

$$\begin{array}{cc}\mathrm{\Omega }^b:=\{(x_1,x_2,\mathrm{},x_{n1})|\hfill & (x_1,x_2,\mathrm{},x_{n1})_{i=1,i<jn2}^{n3}(B_i,B_j),\hfill \\ & B_i\mathrm{\Omega }_{ei}^b,i\{1,2,\mathrm{},n2\}\}.\hfill \end{array}$$

Lemma 2 Suppose $`a(s)=s^n+a_1s^{n1}+\mathrm{}+a_nH^n,b(s)=s^n+b_1s^{n1}+\mathrm{}+b_nH^n,`$ if $`\mathrm{\Omega }^a\mathrm{\Omega }^b\varphi ,`$ take $`(x_1,x_2,\mathrm{},x_{n1})\mathrm{\Omega }^a\mathrm{\Omega }^b,`$ and let $`c(s):=s^{n1}+(x_1\epsilon )s^{n2}+x_2s^{n3}+\mathrm{}+x_{n2}s+(x_{n1}+\epsilon )`$ ($`\epsilon `$ is a sufficiently small positive number), then for $`{\displaystyle \frac{c(s)}{a(s)}}`$ and $`{\displaystyle \frac{c(s)}{b(s)}}`$, we have $`\omega R,\text{Re}[{\displaystyle \frac{c(j\omega )}{a(j\omega )}}]>0`$ and $`\text{Re}[{\displaystyle \frac{c(j\omega )}{b(j\omega )}}]>0.`$

Proof Suppose $`(x_1,x_2,\mathrm{},x_{n1})\mathrm{\Omega }^a\mathrm{\Omega }^b,`$ let $`c(s):=s^{n1}+(x_1\epsilon )s^{n2}+x_2s^{n3}+\mathrm{}+x_{n2}s+(x_{n1}+\epsilon ),\epsilon >0,\epsilon `$ sufficiently small.

$`\omega R,`$ consider

$$\begin{array}{ccc}\text{Re}[\frac{c(j\omega )}{a(j\omega )}]\hfill & =\hfill & \frac{1}{a(j\omega )^2}[c_1\omega ^{2(n1)}+c_2\omega ^{2(n2)}+\mathrm{}+c_{n1}\omega ^2+c_n]\hfill \\ & & +\text{Re}[\frac{\epsilon (j\omega )^{n2}+\epsilon }{a(j\omega )}]\hfill \\ & =\hfill & \frac{1}{a(j\omega )^2}[c_1\omega ^{2(n1)}+c_2\omega ^{2(n2)}+\mathrm{}+c_{n1}\omega ^2+c_n]\hfill \\ & & +\frac{(\epsilon )}{|a(j\omega )|^2}(\omega ^{2(n1)}+\stackrel{~}{c}(\omega ^2)),\hfill \end{array}$$

where $`c_l:=\underset{j=0}{\overset{n}{}}(1)^{l+j}a_jx_{2lj1},l=1,2,\mathrm{},n,`$ $`a_0=1,x_0=1,`$ $`a_i=0`$ if $`i<0`$ or $`i>n,`$ and $`x_i=0`$ if $`i<0`$ or $`i>n1,`$ and $`\stackrel{~}{c}(\omega ^2)`$ is a real polynomial with order not greater than $`2(n2)`$.

In order to prove that $`\omega R,\text{Re}[{\displaystyle \frac{c(j\omega )}{a(j\omega )}}]>0,`$ let $`t=\omega ^2,`$ we only need to prove that, for any $`\epsilon >0,\epsilon `$ sufficiently small, the following polynomial $`f_1(t)`$ satisfies

$$\begin{array}{cc}f_1(t):=\hfill & c_1t^{n1}+c_2t^{n2}+\mathrm{}+c_{n1}t+c_n\hfill \\ & +\epsilon (t^{n1}\stackrel{~}{c}(t))>0,t[0,+\mathrm{}).\hfill \end{array}$$

Since $`(x_1,x_2,\mathrm{},x_{n1})\mathrm{\Omega }^a,`$ by the definition of $`\mathrm{\Omega }^a,`$ it is easy to know that

$$g_1(t):=c_1t^{n1}+c_2t^{n2}+\mathrm{}+c_{n1}t+c_n>0,t(0,+\mathrm{}).$$

Moreover, we obviously have $`f_1(0)>0,`$ and for any $`\epsilon >0`$, when $`t`$ is a sufficiently large or sufficiently small positive number, we have $`f_1(t)>0,`$ namely, there exist $`0<t_1<t_2`$ such that, for all $`\epsilon >0`$, $`t[0,t_1][t_2,+\mathrm{})`$, we have $`f_1(t)>0.`$

Denote

$$M_1=\underset{t[t_1,t_2]}{inf}g_1(t),$$

$$N_1=\underset{t[t_1,t_2]}{sup}|t^{n1}\stackrel{~}{c}(t))|.$$

Then $`M_1>0`$ and $`N_1>0.`$ Choosing $`0<\epsilon <{\displaystyle \frac{M_1}{N_1}},`$ by a direct calculation, we have

$$\begin{array}{cc}f_1(t):=\hfill & c_1t^{n1}+c_2t^{n2}+\mathrm{}+c_{n1}t+c_n\hfill \\ & +\epsilon (t^{n1}\stackrel{~}{c}(t))>0,t[0,+\mathrm{}).\hfill \end{array}$$

Namely,

$$\omega R,\text{Re}[\frac{c(j\omega )}{a(j\omega )}]>0.$$

Similarly, since $`(x_1,x_2,\mathrm{},x_{n1})\mathrm{\Omega }^b,`$ there exist $`0<t_3<t_4`$ such that, for all $`\epsilon >0`$, $`t[0,t_3][t_4,+\mathrm{})`$, we have $`f_2(t)>0,`$where

$$\begin{array}{cc}f_2(t):=\hfill & d_1t^{n1}+d_2t^{n2}+\mathrm{}+d_{n1}t+d_n\hfill \\ & +\epsilon (t^{n1}\stackrel{~}{d}(t)),\hfill \end{array}$$

$$d_l:=\underset{j=0}{\overset{n}{}}(1)^{l+j}b_jx_{2lj1},l=1,2,\mathrm{},n,$$

where $`b_0=1,x_0=1,`$ $`b_i=0`$ if $`i<0`$ or $`i>n,`$ and $`x_i=0`$ if $`i<0`$ or $`i>n1,`$ and $`\stackrel{~}{d}(\omega ^2)`$ is a real polynomial with order not greater than $`2(n2)`$ which is determined by the following equation:

$$\text{Re}[\frac{\epsilon (j\omega )^{n2}+\epsilon }{b(j\omega )}]=\frac{(\epsilon )}{|b(j\omega )|^2}(\omega ^{2(n1)}+\stackrel{~}{d}(\omega ^2)).$$

Denote

$$g_2(t):=d_1t^{n1}+d_2t^{n2}+\mathrm{}+d_{n1}t+d_n,$$

$$M_2=\underset{t[t_3,t_4]}{inf}g_2(t),$$

$$N_2=\underset{t[t_3,t_4]}{sup}|t^{n1}\stackrel{~}{d}(t))|.$$

Then $`M_2>0`$ and $`N_2>0.`$ Choosing $`0<\epsilon <{\displaystyle \frac{M_2}{N_2}},`$ we have

$$\omega R,\text{Re}[\frac{c(j\omega )}{b(j\omega )}]>0.$$

Thus, by choosing $`0<\epsilon <\mathrm{min}\{{\displaystyle \frac{M_1}{N_1}},{\displaystyle \frac{M_2}{N_2}}\},`$ Lemma 2 is proved.

Lemma 3 Suppose $`a(s)=s^n+a_1s^{n1}+\mathrm{}+a_nH^n,b(s)=s^n+b_1s^{n1}+\mathrm{}+b_nH^n,`$ if $`\lambda b(s)+(1\lambda )a(s)H^n,\lambda [0,1],`$ then $`\mathrm{\Omega }^a\mathrm{\Omega }^b\varphi `$

Proof  If $`\lambda [0,1],a_\lambda (s):=\lambda b(s)+(1\lambda )a(s)H^n,`$ by Lemma 1, for any $`\lambda [0,1],`$ $`\mathrm{\Omega }_{ek}^{a_\lambda },k=1,2,\mathrm{},n2,`$ are $`n2`$ ellipses in the first quadrant of the $`R^{n1}`$ space $`(x_1,x_2,\mathrm{},x_{n1})`$.

$`\lambda [0,1],`$ denote

$$\begin{array}{cc}\mathrm{\Omega }^{a_\lambda }:=\{(x_1,x_2,\mathrm{},x_{n1})|\hfill & (x_1,x_2,\mathrm{},x_{n1})_{i=1,i<jn2}^{n3}(A_{\lambda i},A_{\lambda j}),\hfill \\ & A_{\lambda i}\mathrm{\Omega }_{ei}^{a_\lambda },i\{1,2,\mathrm{},n2\}\}\hfill \end{array}$$

Apparently, when $`\lambda `$ changes continuously from $`0`$ to $`1`$, $`\mathrm{\Omega }^{a_\lambda }`$ will change continuously from $`\mathrm{\Omega }^a`$ to $`\mathrm{\Omega }^b,`$ and $`\mathrm{\Omega }_{ek}^{a_\lambda }`$ will change continuously from $`\mathrm{\Omega }_{ek}^a`$ to $`\mathrm{\Omega }_{ek}^b,k=1,2,\mathrm{},n2.`$

Now assume $`\mathrm{\Omega }^a\mathrm{\Omega }^b=\varphi ,`$ by the definitions of $`\mathrm{\Omega }^a`$ and $`\mathrm{\Omega }^b`$ and Lemma 1, $`u_1>0,`$ $`u_2>0,u_1a_1,u_1b_1,`$ and $`\stackrel{~}{k}\{1,2,\mathrm{},n2\},`$ such that the following hyperplane $`L`$ in the $`R^{n1}`$ space $`(x_1,x_2,\mathrm{},x_{n1})`$

$$L:\frac{x_1}{u_1}+\frac{x_2}{u_2}+\mathrm{}+\frac{x_{n1}}{u_{n1}}=1$$

separates $`\mathrm{\Omega }^a`$ and $`\mathrm{\Omega }^b,`$ meanwhile, $`L`$ is tangent with $`\mathrm{\Omega }_{e1}^a,\mathrm{\Omega }_{e2}^a,\mathrm{},\mathrm{\Omega }_{e(n2)}^a`$ and $`\mathrm{\Omega }_{e\stackrel{~}{k}}^b`$ simultaneously (or tangent with $`\mathrm{\Omega }_{e1}^b,\mathrm{\Omega }_{e2}^b,\mathrm{},\mathrm{\Omega }_{e(n2)}^b`$ and $`\mathrm{\Omega }_{e\stackrel{~}{k}}^a`$ simultaneously).

Without loss of generality, suppose that $`L`$ is tangent with $`\mathrm{\Omega }_{e1}^a,\mathrm{\Omega }_{e2}^a,\mathrm{},\mathrm{\Omega }_{e(n2)}^a`$ and $`\mathrm{\Omega }_{e\stackrel{~}{k}}^b`$ simultaneously.

In what follows, the notation $`[x]`$ stands for the largest integer that is smaller than or equal to the real number $`x`$, and $`y_z`$ stands for the remainder of the nonnegative integer $`y`$ divided by the positive integer $`z`$ 333 For example, $`[1.5]=1,[0.5]=0,[1.5]=2`$, and $`0_2=0,1_2=1,11_3=2`$..

Since $`L`$ is tangent with $`\mathrm{\Omega }_{e1}^a,\mathrm{\Omega }_{e2}^a,\mathrm{},\mathrm{\Omega }_{e(n2)}^a`$ and $`\mathrm{\Omega }_{e\stackrel{~}{k}}^b`$ simultaneously, note that $`a(s)`$ is Hurwitz stable and $`u_1>0,u_1a_1,u_2>0,`$ using mathematical induction, by a lengthy calculation, we know that the necessary and sufficient condition for $`L`$ being tangent with $`\mathrm{\Omega }_{e1}^a,\mathrm{\Omega }_{e2}^a,\mathrm{},\mathrm{\Omega }_{e(n2)}^a`$ simultaneously is 444 When $`n=3,`$ we have (see  for details):
$$u_1u_2a_1u_2a_2u_1+a_3=0.$$
When $`n=4,`$ we have (see  - for details):
$$\{\begin{array}{c}u_1u_2^2a_1u_2^2a_2u_1u_2+a_3u_2+a_4u_1=0,\hfill \\ u_3=u_1u_2.\hfill \end{array}$$
When $`n=5,`$ we have (see  for details):
$$\{\begin{array}{c}u_1u_2^2a_1u_2^2a_2u_1u_2+a_3u_2+a_4u_1a_5=0,\hfill \\ u_3=u_1u_2,u_4=u_2^2.\hfill \end{array}$$
When $`n=6,`$ we have (see  for details):
$$\{\begin{array}{c}u_1u_2^3a_1u_2^3a_2u_1u_2^2+a_3u_2^2+a_4u_1u_2a_5u_2a_6u_1=0,\hfill \\ u_3=u_1u_2,u_4=u_2^2,u_5=u_1u_2^2.\hfill \end{array}$$

$$\underset{i=0}{\overset{n}{}}(1)^{\left[{\displaystyle \frac{i+1}{2}}\right]}a_iu_1^{i+1_2}u_2^{\left[{\displaystyle \frac{n}{2}}\right]\left[{\displaystyle \frac{i}{2}}\right]}=0$$
(1)
and

$$u_j=(1)^{\left[{\displaystyle \frac{j1}{2}}\right]}u_1^{j_2}u_2^{\left[{\displaystyle \frac{j}{2}}\right]},j=3,4,\mathrm{},n1,$$
(2)
where $`a_0=1.`$

Since $`u_j=(1)^{\left[{\displaystyle \frac{j1}{2}}\right]}u_1^{j_2}u_2^{\left[{\displaystyle \frac{j}{2}}\right]},j=3,4,\mathrm{},n1,L`$ is tangent with $`\mathrm{\Omega }_{e\stackrel{~}{k}}^b`$, by a direct calculation, we have

$$\underset{i=0}{\overset{n}{}}(1)^{\left[{\displaystyle \frac{i+1}{2}}\right]}b_iu_1^{i+1_2}u_2^{\left[{\displaystyle \frac{n}{2}}\right]\left[{\displaystyle \frac{i}{2}}\right]}=0$$
(3)
where $`b_0=1.`$

From $`(\text{1})`$ and $`(\text{3}),`$ we obviously have $`\lambda [0,1],`$

$$\underset{i=0}{\overset{n}{}}(1)^{\left[{\displaystyle \frac{i+1}{2}}\right]}a_{\lambda i}u_1^{i+1_2}u_2^{\left[{\displaystyle \frac{n}{2}}\right]\left[{\displaystyle \frac{i}{2}}\right]}=0$$
(4)
where $`a_{\lambda i}:=a_i+\lambda (b_ia_i),a_0=1,b_0=1,i=0,1,2,\mathrm{},n.`$ $`(\text{4})`$ and $`(\text{2})`$ show that $`L`$ is also tangent with $`\mathrm{\Omega }_{e\stackrel{~}{k}}^{a_\lambda }(\lambda [0,1])`$, but $`L`$ separates $`\mathrm{\Omega }_{e\stackrel{~}{k}}^a`$ and $`\mathrm{\Omega }_{e\stackrel{~}{k}}^b,`$ and when $`\lambda `$ changes continuously from $`0`$ to $`1`$, $`\mathrm{\Omega }_{e\stackrel{~}{k}}^{a_\lambda }`$ will change continuously from $`\mathrm{\Omega }_{e\stackrel{~}{k}}^a`$ to $`\mathrm{\Omega }_{e\stackrel{~}{k}}^b,`$ which is obviously impossible. This completes the proof.

Lemma 4 Suppose $`a(s)=s^n+a_1s^{n1}+\mathrm{}+a_nH^n,b(s)=s^n+b_1s^{n1}+\mathrm{}+b_nH^n,c(s)=s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1},`$ if $`\omega R,\text{Re}[{\displaystyle \frac{c(j\omega )}{a(j\omega )}}]>0`$ and $`\text{Re}[{\displaystyle \frac{c(j\omega )}{b(j\omega )}}]>0,`$ take

$`\stackrel{}{c}(s):=c(s)+\delta h(s),\delta >0,\delta `$ sufficiently small,

(where $`h(s)`$ is an arbitrarily given monic $`n`$-$`th`$ order polynomial), then $`{\displaystyle \frac{\stackrel{}{c}(s)}{a(s)}}`$ and $`{\displaystyle \frac{\stackrel{}{c}(s)}{b(s)}}`$ are both strictly positive real.

Proof  Obviously, $`(\stackrel{~}{c})=(a)=n,`$ namely, $`\stackrel{~}{c}(s)`$ and $`a(s)`$ have the same order $`n`$. Since $`a(s)H^n,`$ there exists $`\omega _1>0`$ such that, for all $`\omega \omega _1,`$ we have Re$`({\displaystyle \frac{\stackrel{~}{c}(j\omega )}{a(j\omega )}})>0.`$

Denote

$$M_3=\underset{\omega \omega _1}{inf}\text{Re}(\frac{c(j\omega )}{a(j\omega )}),N_3=\underset{\omega \omega _1}{sup}\text{Re}(\frac{h(j\omega )}{a(j\omega )}).$$

Then $`M_3>0`$ and $`N_3>0.`$ Choosing $`0<\epsilon <{\displaystyle \frac{M_3}{N_3}},`$ it can be directly verified that

$$\text{Re}(\frac{\stackrel{~}{c}(j\omega )}{a(j\omega )})>0,\omega R.$$

Similarly, $`(\stackrel{~}{c})=(b)=n,`$ and there exists $`\omega _2>0`$ such that, for all $`\omega \omega _2,`$ we have Re$`({\displaystyle \frac{\stackrel{~}{c}(j\omega )}{b(j\omega )}})>0.`$

Denote

$$M_4=\underset{\omega \omega _2}{inf}\text{Re}(\frac{c(j\omega )}{b(j\omega )}),N_4=\underset{\omega \omega _2}{sup}\text{Re}(\frac{h(j\omega )}{b(j\omega )}).$$

Then $`M_4>0`$ and $`N_4>0.`$ Choosing $`0<\epsilon <{\displaystyle \frac{M_4}{N_4}},`$ it can be directly verified that

$$\text{Re}(\frac{\stackrel{~}{c}(j\omega )}{b(j\omega )})>0,\omega R.$$

Thus, by choosing $`0<\epsilon <\mathrm{min}\{{\displaystyle \frac{M_3}{N_3}},{\displaystyle \frac{M_4}{N_4}}\},`$ Lemma 4 is proved.

The sufficiency of Theorem 1 is now proved by simply combining Lemmas 1-4.

Remark 1 From the proof of Theorem 1, we can see that this paper not only proves the existence, but also provides a design method. In fact, based on the main idea of this paper, we have developed a geometric algorithm with order reduction for robust SPR synthesis which is very efficient for high order polynomial segments .

Remark 2 The method provided in this paper is constructive, and is insightful and helpful in solving more general robust SPR synthesis problems for polynomial polytopes, multilinear families, etc..

Remark 3 Our main results in this paper can also be extended to discrete-time case. In fact, by applying the bilinear transformation, we can transform the unit circle into the left half plane. Hence, Theorem 1 can be generalized to discrete-time case. Moreover, in discrete-time case, the order of the polynomial obtained by our method is bounded by the order of the given polynomial segment .

Remark 4 If $`{\displaystyle \frac{c(s)}{a(s)}}`$ and $`{\displaystyle \frac{c(s)}{b(s)}}`$ are both SPR, it is easy to know that $`\lambda [0,1],{\displaystyle \frac{c(s)}{\lambda a(s)+(1\lambda )b(s)}}`$ is also SPR.

Remark 5 The stability of polynomial segment can be checked by many efficient methods, e.g., eigenvalue method, root locus method, value set method, etc. .

## 3 Conclusions

We have constructively proved that, for any two $`n`$-$`th`$ order polynomials $`a(s)`$ and $`b(s),`$ the Hurwitz stability of their convex combination is necessary and sufficient for the existence of a polynomial $`c(s)`$ such that $`c(s)/a(s)`$ and $`c(s)/b(s)`$ are both strictly positive real. By using similar method, we can also constructively prove the existence of SPR synthesis for low order ($`n4`$) interval polynomials. Namely, when $`n4,`$ the Hurwitz stability of the four Kharitonov vertex polynomials is necessary and sufficient for the existence of a fixed polynomial such that the ratio of this polynomial to any polynomial in the interval polynomial set is SPR invariant . The SPR synthesis problem for high order interval polynomials is currently under investigation.


text:
# Approximations to the QED Fermion Green’s Function in a Constant External Field

## I Introduction

An exact representation for the causal QED fermion Green’s function, $`G_c(x,y|A)`$, in an arbitrary external field, was derived by Fried et al. in such a way as to obtain an exact representation which naturally allows for non-perturbative approximations. The rather intimidating exact representation is

$`G_c(x,y|A)=`$ $`i`$ $`{\displaystyle _0^{\mathrm{}}}𝑑se^{ism^2}{\displaystyle \frac{d^4p}{(2\pi )^4}e^{ip(xy)}\frac{d^4zd^4P}{(2\pi )^4}e^{iP(zy)+i\frac{s}{4}P^2}}`$ (1)
$`\times `$ $`{\displaystyle \underset{N=1}{\overset{\mathrm{}}{}}^{}}{\displaystyle \frac{(i)^2}{(2\pi )^4}}{\displaystyle d^4P_Nd^4Q_Ne^{\frac{i}{2}(P_N^2+Q_N^2)}}`$ (2)
$`\times `$ $`e^{i_0^s𝑑s^{}[p\mathrm{\Omega }(s^{})]^2}e^{g_0^s𝑑s^{}\frac{}{z_\mu }A_\mu \left(\zeta (s^{})2_0^s^{}\mathrm{\Omega }\right)}`$ (3)
$`\times `$ $`\{mi\gamma [p\mathrm{\Omega }(s)]\}\left(e^{g_0^s𝑑s^{}\sigma F\left(\zeta (s^{})2_0^s^{}\mathrm{\Omega }\right)}\right)_+,`$ (4)
where $`\mathrm{\Omega }(s^{})`$ is the solution to the “map”

$$\mathrm{\Omega }(s^{})=gA\left(\zeta (s^{})2_0^s^{}\mathrm{\Omega }\right),$$
(5)
and $`\zeta (s^{})`$ is given by

$$\zeta (s^{})=z+s^{}(2p+P)\frac{2\sqrt{s}}{\pi }\underset{N=1}{\overset{\mathrm{}}{}}{}_{}{}^{}\frac{1}{N}\left[P_N\mathrm{cos}\left(\frac{N\pi s^{}}{s}\right)+Q_N\mathrm{sin}\left(\frac{N\pi s^{}}{s}\right)\right].$$
(6)
$`_{N=1}^{\mathrm{}}{}_{}{}^{}`$ and $`_{N=1}^{\mathrm{}}{}_{}{}^{}`$ represent the product and sum respectively over all odd natural numbers. It is easy to see that we may approximate the Green’s function non-perturbatively by retaining a finite number of integrations, in particular, for the $`n`$th approximation $`G_c^{(n)}(x,y|A)`$, where $`n=0,1,2,\mathrm{}`$, we retain $`2n+3`$ integrations. This means that for $`n0`$, the odd number $`N`$ takes the values $`N=1,3,\mathrm{},2n1`$, while for the zeroth approximation all $`N`$-dependence is neglected. The exact result is of course recovered in the limit $`n\mathrm{}`$, that is, $`G_c(x,y|A)=lim_n\mathrm{}G_c^{(n)}(x,y|A)`$.

It is well known that there is an exact expression, first obtained by Schwinger , for the fermion Green’s function in the simple case of a constant, but otherwise arbitrary, external field, which (in an arbitrary gauge) reads

$`G_c(x,y|A)=`$ $`\mathrm{\Phi }(x,y|A){\displaystyle \frac{1}{(4\pi )^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{ds}{s^2}}e^{ism^2}e^{gs\sigma F}\left(det{\displaystyle \frac{\mathrm{sinh}gFs}{gFs}}\right)^{\frac{1}{2}}`$ (7)
$`\times `$ $`e^{\frac{i}{4}(xy)gF\mathrm{coth}gFs(xy)}\left[m{\displaystyle \frac{i}{2}}\gamma (gF\mathrm{coth}gFs+gF)(xy)\right],`$ (8)
where the holonomy factor $`\mathrm{\Phi }(x,y|A)`$ is

$$\mathrm{\Phi }(x,y|A)=e^{ig_y^x𝑑\xi _\mu [A_\mu (\xi )+\frac{1}{2}F(\xi y)]},$$
(9)
and carries the complete gauge-dependence of the Green’s function. It is convenient to employ matrix notation, in which we regard the field strength tensor as a constant, antisymmetric $`4\times 4`$ matrix $`F`$.

It is obvious that Schwinger’s result must somehow be contained as a special case of the exact representation (4)-(6). Furthermore, the Fradkin representation ,, from which the above representation was derived, almost trivially yields Schwinger’s result. This leads us to expect that the latter may be extracted, analytically, from the above representation, and here we demonstrate that this is indeed the case.

We proceed by evaluating that $`n`$th non-perturbative approximation to the Green’s function in a constant field as given by (4)-(6) Schwinger’s result is then recovered in the limit as the order of this approximation approaches infinity. That we may carry out this programme relies on two things: for the case of constant $`F`$ the ordered exponential of (4) becomes an ordinary exponential, and the $`2n+3`$ integrations we must perform are all Gaussian. Indeed a large part of the evaluation of the $`n`$th approximation is an extension of our matrix notation to account for this latter fact. $`2n+3`$ Gaussian integrals may be expressed as one Gaussian integral over a $`2n+3`$-dimensional space. The integration is then, with the appropriate notation, essentially trivial. In Section II we make the change of notation, perform the integration, and express the $`n`$th approximation in terms of the two resulting, order $`2n+3`$, determinants. In Section III we evaluate the determinants and take the limit $`n\mathrm{}`$ to recover Schwinger’s result.

We emphasise that the result when stopped at a finite $`n`$ (except $`n=0`$) is essentially non-perturbative, in that it is *not* a polynomial in $`g`$. The result of (4)-(6) is thus, as $`n`$ is increased, a systematically improving, but non-perturbative, method for calculating the Green’s function in an external field.

## II The Integration

We must first choose a gauge to work in. It shall be convenient to work exclusively in the Schwinger-Fock ($`SF`$) gauge

$$A_\mu ^{SF}(z)=\frac{1}{2}F_{\mu \nu }(zy)_\nu ,$$
(10)
the initial motivation for which is that we may forget about the holonomy factor, which reduces to 1. An immediate conseqence is that $`\frac{}{z_\mu }A_\mu ^{SF}(z)=0`$, due to the antisymmetry of $`F`$; the factor in (4) containing this term in the exponent also reduces to 1. With the above simplifications, the $`n`$th approximation to the fermion Green’s function in a constant field, in the $`SF`$ gauge, as given by (4), is thus

$`G_c^{(n)}(x,y|A^{SF})=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑se^{ism^2}e^{gs\sigma F}{\displaystyle \frac{i}{(2\pi )^8}}{\displaystyle \frac{(i)^{2n}}{(2\pi )^{4n}}}{\displaystyle 𝑑p𝑑z𝑑P\underset{N=1}{\overset{2n1}{}}^{}dP_NdQ_N}`$ (11)
$`\times `$ $`e^{ip(xy)+iP(zy)+i\frac{s}{4}P^2+\underset{N=1}{\overset{2n1}{}}{}_{}{}^{}\frac{i}{2}(P_N^2+Q_N^2)}`$ (12)
$`\times `$ $`e^{i_0^s𝑑s^{}[p\mathrm{\Omega }^{(n)}(s^{})]^2}\{mi\gamma [p\mathrm{\Omega }^{(n)}(s)]\},`$ (13)
where $`\mathrm{\Omega }^{(n)}(s^{})`$ is the solution to (5), (6), but with the sum in (6) terminating at $`N=2n1`$, and in the spirit of using matrix notation we have used $`dp`$ instead of $`d^4p`$ and so on. In the $`SF`$ gauge, (5), (6) yield a simple integral equation for $`\mathrm{\Omega }^{(n)}(s^{})`$:

$`\mathrm{\Omega }^{(n)}(s^{})gF{\displaystyle _0^s^{}}\mathrm{\Omega }^{(n)}`$ (14)
$`=`$ $`{\displaystyle \frac{1}{2}}gF\left\{zy+s^{}(2p+P){\displaystyle \frac{2\sqrt{s}}{\pi }}{\displaystyle \underset{N=1}{\overset{2n1}{}}}{}_{}{}^{}{\displaystyle \frac{1}{N}}\left[P_N\mathrm{cos}\left({\displaystyle \frac{N\pi s^{}}{s}}\right)+Q_N\mathrm{sin}\left({\displaystyle \frac{N\pi s^{}}{s}}\right)\right]\right\}.`$ (15)
The equivalent differential equation plus boundary condition may be solved with only elementary integrals. It is the combination $`p\mathrm{\Omega }^{(n)}(s^{})`$ which appears in (13), we find

$`p\mathrm{\Omega }^{(n)}(s^{})=`$ $`e`$ $`{}_{}{}^{gFs^{}}p+{\displaystyle \frac{1}{2}}gFe^{gFs^{}}(zy)+{\displaystyle \frac{1}{2}}(e^{gFs^{}}1)P`$ (16)
$``$ $`{\displaystyle \frac{1}{\sqrt{s}}}{\displaystyle \underset{N=1}{\overset{2n1}{}}}{}_{}{}^{}{\displaystyle \frac{1}{1+\mathrm{\Lambda }_N^2}}\left[\mathrm{sin}\left({\displaystyle \frac{N\pi s^{}}{s}}\right)+\mathrm{\Lambda }_N\mathrm{cos}\left({\displaystyle \frac{N\pi s^{}}{s}}\right)+{\displaystyle \frac{1}{\mathrm{\Lambda }_N}}e^{gFs^{}}\right]P_N`$ (17)
$``$ $`{\displaystyle \frac{1}{\sqrt{s}}}{\displaystyle \underset{N=1}{\overset{2n1}{}}}{}_{}{}^{}{\displaystyle \frac{1}{1+\mathrm{\Lambda }_N^2}}\left[\mathrm{cos}\left({\displaystyle \frac{N\pi s^{}}{s}}\right)+\mathrm{\Lambda }_N\mathrm{sin}\left({\displaystyle \frac{N\pi s^{}}{s}}\right)+e^{gFs^{}}\right]Q_N,`$ (18)
where $`\mathrm{\Lambda }_N=\frac{N\pi }{gFs}`$. Our choice of gauge allows us to make the change of variable $`zyz`$, after which $`p\mathrm{\Omega }^{(n)}(s^{})`$ is independent of $`y`$, so that the only $`x`$ and $`y`$ dependence in the exponent of (13) appears in the term $`ip(xy)`$.

If we imagine substituting (18) into (13), we recognise that all of the terms in the exponent of the integrand, save the term $`ip(xy)`$, are able to be expressed in the form $`\frac{i}{2}X_i^TA_{ij}^{(n)}X_j`$, where the $`X_i`$ are the $`2n+3`$ 4-vector variables we must integrate over, the $`A_{ij}^{(n)}`$ are some matrix functions of $`F`$, and the $`\frac{i}{2}`$ is a convenient normalisation factor. This suggests that we extend our matrix notation, and write this part of the exponent of (13) as $`\frac{i}{2}X^TA^{(n)}X`$, where $`X`$ is a $`(2n+3)\times 1`$ column vector of 4-vector variables, and $`A^{(n)}`$ is a $`(2n+3)\times (2n+3)`$ symmetric matrix with matrix elements $`A_{ij}^{(n)}`$. Now $`dpdzdP\underset{N=1}{\overset{2n1}{}}{}_{}{}^{}dP_NdQ_N=d^{2n+3}X`$. Let us define the column vector $`X`$ such that

$$X^T=[pzPP_1Q_1\mathrm{}P_NQ_N\mathrm{}P_{2n1}Q_{2n1}].$$
(19)
The term $`ip(xy)`$ is included by introducing

$$B^{(n)T}=[(xy)\mathrm{\hspace{0.33em}\hspace{0.33em}0\hspace{0.33em}\hspace{0.33em}0}\mathrm{}\mathrm{\hspace{0.33em}\hspace{0.33em}0}],$$
(20)
whence $`e^{ip(xy)}=e^{iB^{(n)T}X}`$. A construction for the matrix $`A^{(n)}`$ is obtained by defining

$$p\mathrm{\Omega }^{(n)}(s^{})=C^{(n)T}(s^{})X,$$
(21)
so that $`e^{i_0^s𝑑s^{}[p\mathrm{\Omega }^{(n)}(s^{})]^2}=e^{iX^T_0^s𝑑s^{}C^{(n)}(s^{})C^{(n)T}(s^{})X}`$. The elements of the row vector $`C^{(n)T}(s^{})`$ are read straight from (18). The transpose of the row vector is

$$C^{(n)}(s^{})=\left[\begin{array}{c}e^{gFs^{}}\\ \frac{1}{2}gFe^{gFs^{}}\\ \frac{1}{2}(e^{gFs^{}}1)\\ \frac{1}{\sqrt{s}}\frac{1}{1+\mathrm{\Lambda }_1^2}\left[\mathrm{sin}\left(\frac{\pi s^{}}{s}\right)+\mathrm{\Lambda }_1\mathrm{cos}\left(\frac{\pi s^{}}{s}\right)+\frac{1}{\mathrm{\Lambda }_1}e^{gFs^{}}\right]\\ \frac{1}{\sqrt{s}}\frac{1}{1+\mathrm{\Lambda }_1^2}\left[\mathrm{cos}\left(\frac{\pi s^{}}{s}\right)+\mathrm{\Lambda }_1\mathrm{sin}\left(\frac{\pi s^{}}{s}\right)e^{gFs^{}}\right]\\ \mathrm{}\\ \frac{1}{\sqrt{s}}\frac{1}{1+\mathrm{\Lambda }_N^2}\left[\mathrm{sin}\left(\frac{N\pi s^{}}{s}\right)+\mathrm{\Lambda }_N\mathrm{cos}\left(\frac{N\pi s^{}}{s}\right)+\frac{1}{\mathrm{\Lambda }_N}e^{gFs^{}}\right]\\ \frac{1}{\sqrt{s}}\frac{1}{1+\mathrm{\Lambda }_N^2}\left[\mathrm{cos}\left(\frac{N\pi s^{}}{s}\right)+\mathrm{\Lambda }_N\mathrm{sin}\left(\frac{N\pi s^{}}{s}\right)e^{gFs^{}}\right]\\ \mathrm{}\end{array}\right],$$
(22)
where we have used $`F^T=F`$ $`(\mathrm{\Lambda }_N^T=\mathrm{\Lambda }_N)`$, and omitted the last two $`N=2n1`$ elements for brevity. We also write $`e^{iPz+i\frac{s}{4}P^2+\underset{N=1}{\overset{2n1}{}}{}_{}{}^{}\frac{i}{2}(P_N^2+Q_N^2)}=e^{\frac{i}{2}X^TD^{(n)}X}`$, which defines the matrix $`D^{(n)}`$:

$$D^{(n)}=\left[\begin{array}{cccccccc}0& 0& 0& 0& 0& \mathrm{}& 0& 0\\ 0& 0& 1& 0& 0& \mathrm{}& 0& 0\\ 0& 1& \frac{s}{2}& 0& 0& \mathrm{}& 0& 0\\ 0& 0& 0& 1& 0& \mathrm{}& 0& 0\\ 0& 0& 0& 0& 1& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}& 1& 0\\ 0& 0& 0& 0& 0& \mathrm{}& 0& 1\end{array}\right].$$
(23)
Our change of notation complete, (13) and (18) become

$`G_c^{(n)}(x,y|A^{SF})=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑se^{ism^2}e^{gs\sigma F}{\displaystyle \frac{i}{(2\pi )^8}}{\displaystyle \frac{(i)^{2n}}{(2\pi )^{4n}}}`$ (24)
$`\times `$ $`{\displaystyle d^{2n+3}Xe^{\frac{i}{2}X^TA^{(n)}X+iB^{(n)T}X}[mi\gamma C^{(n)T}(s)X]},`$ (25)
where

$$A^{(n)}=2_0^s𝑑s^{}C^{(n)}(s^{})C^{(n)T}(s^{})+D^{(n)}.$$
(26)
All that we have done is to re-express, in the usual way, the product of $`2n+3`$ Gaussian integrals as one $`(2n+3)`$-dimensional Gaussian integral. After a change of variable $`XXA^{(n)1}B^{(n)}`$, the integration is trivial, we obtain

$`G_c^{(n)}(x,y|A^{SF})=`$ $`{\displaystyle \frac{1}{4\pi ^2}}`$ $`{\displaystyle _0^{\mathrm{}}}𝑑se^{ism^2}e^{gs\sigma F}(\text{Det}A^{(n)})^{\frac{1}{2}}e^{\frac{i}{2}B^{(n)T}A^{(n)1}B^{(n)}}`$ (27)
$`\times `$ $`[m+i\gamma C^{(n)T}(s)A^{(n)1}B^{(n)}],`$ (28)
where $`\text{Det}A^{(n)}`$ is the determinant of $`A^{(n)}`$.

The matrix $`A^{(n)}`$is a $`4(2n+3)\times 4(2n+3)`$ matrix which is naturally partitioned into the $`(2n+3)\times (2n+3)`$ $`4\times 4`$ matrices $`A_{ij}^{(n)}=2_0^s𝑑s^{}C_i^{(n)}(s^{})C_j^{(n)T}(s^{})+D_{ij}^{(n)}`$. The $`A_{ij}^{(n)}`$, as (matrix) functions of $`F`$ only, commute, and are thus referred to as the “elements” of $`A^{(n)}`$. Since $`A^{(n)}`$ is symmetric, we have the relations $`A_{ii}^{(n)}=A_{ii}^{(n)T}`$ and $`A_{ji}^{(n)}=A_{ij}^{(n)T}`$. If we form a determinant using the elements $`A_{ij}^{(n)}`$, the result will be a $`4\times 4`$ matrix, also some function of $`F`$, which we denote $`detA^{(n)}`$. It is easy to convince oneself that the determinant of $`A^{(n)}`$, which may be partitioned in this way, is the determinant of the matrix we have called $`detA^{(n)}`$, that is, $`\text{Det}A^{(n)}=det(detA^{(n)})`$.

Recalling that $`B_i^{(n)}=(xy)\delta _{i1}`$, and with the notation discussed above, we may write

$$C^{(n)T}(s)A^{(n)1}B^{(n)}=\frac{detA_1^{(n)}}{detA^{(n)}}(xy),$$
(29)
where $`A_1^{(n)}`$ is the matrix obtained by replacing the first row of $`A^{(n)}`$ by $`C^{(n)T}(s)`$. Similarly,

$$B^{(n)T}A^{(n)1}B^{(n)}=(xy)\frac{det\stackrel{~}{A}^{(n)}}{detA^{(n)}}(xy),$$
(30)
where $`det\stackrel{~}{A}^{(n)}`$ is the $`(1,1)`$ cofactor of $`A^{(n)}`$, that is, $`\stackrel{~}{A}^{(n)}`$ is the matrix obtained by deleting the first row and the first column of $`A^{(n)}`$. In fact it will be shown in the next section that $`\frac{det\stackrel{~}{A}^{(n)}}{detA^{(n)}}`$ is the part of $`\frac{detA_1^{(n)}}{detA^{(n)}}`$ symmetric with respect to the interchange of space-time indices, for all $`n`$, and denoted by a superscript $`S`$. The $`n`$th approximation to the Green’s function, in terms of the two determinants $`detA^{(n)}`$ and $`detA_1^{(n)}`$, is then

$`G_c^{(n)}(x,y|A^{SF})=`$ $`{\displaystyle \frac{1}{4\pi ^2}}`$ $`{\displaystyle _0^{\mathrm{}}}𝑑se^{ism^2}e^{gs\sigma F}(\text{Det}A^{(n)})^{\frac{1}{2}}e^{\frac{i}{2}(xy)\left(\frac{detA_1^{(n)}}{detA^{(n)}}\right)^S(xy)}`$ (31)
$`\times `$ $`\left[m+i\gamma {\displaystyle \frac{detA_1^{(n)}}{detA^{(n)}}}(xy)\right].`$ (32)
Comparing (32) with Schwinger’s result (8), and letting $`detA=lim_n\mathrm{}detA^{(n)}`$ and $`detA_1=lim_n\mathrm{}detA_1^{(n)}`$, it is necessary that

$$detA=2s\frac{\mathrm{sinh}gFs}{gFs},$$
(33)
so that $`(\text{Det}A)^{\frac{1}{2}}=\frac{1}{4s^2}(det\frac{\mathrm{sinh}gFs}{gFs})^{\frac{1}{2}}`$, and

$$detA_1=e^{gFs},$$
(34)
so that the quotient $`\frac{detA_1}{detA}=\frac{1}{2}(gF\mathrm{coth}gFs+gF)`$. Note that the symmetric part of $`\frac{detA_1}{detA}`$ is $`(\frac{detA_1}{detA})^S=\frac{1}{2}gF\mathrm{coth}gFs`$, we mentioned above that this relation holds for all values of $`n`$. Thus $`detA^{(n)}`$ and $`detA_1^{(n)}`$ are yet to be determined, $`n`$th, non-perturbative approximations to (33) and (34) respectively. The quotient $`\frac{detA_1^{(n)}}{detA^{(n)}}`$ provides a non-perturbative approximation to $`\frac{1}{2}(gF\mathrm{coth}gFs+gF)`$, the symmetric part of the former a non-perturbative approximation to the symmetric part of the latter. In the next section we calculate exact expressions for these non-perturbative approximations, and show that we can obtain (33) and (34), and thus Schwinger’s result, in the limit $`n\mathrm{}`$.

## III The Determinants

An obvious but important fact is that we only need to find the determinants of the matrices $`A^{(n)}`$ and $`A_1^{(n)}`$, not the matrices themselves. This means that we can simplify $`C^{(n)}(s^{})`$, its transpose, and $`D^{(n)}`$, with any row and column operations which do not alter the determinant, before using (34) to find $`A^{(n)}`$. We now find the reduced form of $`A^{(n)}`$, from this it will be easy to obtain the reduced form of $`A_1^{(n)}`$ and of $`\stackrel{~}{A}^{(n)}`$. We perform the following sets of row operations on the column vector $`C^{(n)}(s^{})`$: use the first element to remove those terms proportional to $`e^{gFs^{}}`$ from all other elements, noting that the second element requires the row operation row2 $``$ row$`{}_{2}{}^{}+\frac{1}{2}gF`$row1; then rowN+3 $``$ row$`{}_{N+3}{}^{}\mathrm{\Lambda }_N`$rowN+4; then row$`{}_{N+4}{}^{}`$ row$`{}_{N+4}{}^{}+\frac{\mathrm{\Lambda }_N\frac{1}{\mathrm{\Lambda }_N}}{1+\mathrm{\Lambda }_N^2}`$rowN+3; the last two sets for all rows $`N=1,3,\mathrm{},2n1`$. $`C^{(n)}(s^{})`$ becomes

$$C^{(n)}(s^{})\left[\begin{array}{c}e^{gFs^{}}\\ 0\\ \frac{1}{2}\\ \frac{1}{\sqrt{s}}\mathrm{sin}\left(\frac{\pi s^{}}{s}\right)\\ \frac{1}{\sqrt{s}}\frac{1}{1+\mathrm{\Lambda }_1^2}\left[\mathrm{cos}\left(\frac{\pi s^{}}{s}\right)+\frac{1}{\mathrm{\Lambda }_1}\mathrm{sin}\left(\frac{\pi s^{}}{s}\right)\right]\\ \mathrm{}\\ \frac{1}{\sqrt{s}}\mathrm{sin}\left(\frac{N\pi s^{}}{s}\right)\\ \frac{1}{\sqrt{s}}\frac{1}{1+\mathrm{\Lambda }_N^2}\left[\mathrm{cos}\left(\frac{N\pi s^{}}{s}\right)+\frac{1}{\mathrm{\Lambda }_N}\mathrm{sin}\left(\frac{N\pi s^{}}{s}\right)\right]\\ \mathrm{}\end{array}\right].$$
(35)
The row operations are performed on $`C^{(n)}(s^{})`$ and $`D^{(n)}`$. To keep things symmetric we perform the transposed operations on $`C^{(n)T}(s^{})`$ and $`D^{(n)}`$. The integrals required by (34) are elementary and we easily obtain the reduced form of the matrix $`A^{(n)}`$:

$$A^{(n)}\left[\begin{array}{ccccccccc}2s& 0& \frac{1}{gF}(e^{gFs}1)& x_1^T& 0& \mathrm{}& x_N^T& 0& \mathrm{}\\ 0& 0& 1& 0& 0& \mathrm{}& 0& 0& \mathrm{}\\ \frac{1}{gF}(e^{gFs}1)& 1& 0& 0& 0& \mathrm{}& 0& 0& \mathrm{}\\ x_1& 0& 0& \mathrm{\Lambda }_1^2& \frac{3\mathrm{\Lambda }_1}{1+\mathrm{\Lambda }_1^2}& \mathrm{}& 0& 0& \mathrm{}\\ 0& 0& 0& \frac{3\mathrm{\Lambda }_1}{1+\mathrm{\Lambda }_1^2}& \frac{5\mathrm{\Lambda }_1^2}{(1+\mathrm{\Lambda }_1^2)^2}& \mathrm{}& 0& 0& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \\ x_N& 0& 0& 0& 0& \mathrm{}& \mathrm{\Lambda }_N^2& \frac{3\mathrm{\Lambda }_N}{1+\mathrm{\Lambda }_N^2}& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}& \frac{3\mathrm{\Lambda }_N}{1+\mathrm{\Lambda }_N^2}& \frac{5\mathrm{\Lambda }_N^2}{(1+\mathrm{\Lambda }_N^2)^2}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right],$$
(36)
where $`x_N=\frac{2}{\sqrt{s}}\frac{\mathrm{\Lambda }_N}{1+\mathrm{\Lambda }_N^2}\frac{1}{gF}(1+e^{gFs})`$ and we have further used the $`(2,3)`$ and $`(3,2)`$ elements of the reduced $`A^{(n)}`$ to eliminate those elements to the right of the former and below the latter. After trivially expanding the determinant along the second row and down the second column we are left with the determinant of a $`(2n+1)\times (2n+1)`$ bordered matrix. The first set of row and column operations and the orthogonality of the sine and cosine functions have ensured the matrix is block diagonal, the other sets of operations that every second element of the border is zero. Expanding along the first row and down the first column we obtain the following expression for the $`n`$th non-perturbative approximation to (33):

$$detA^{(n)}=2s\left(\underset{N=1}{\overset{2n1}{}}{}_{}{}^{}\alpha _{N}^{}\right)\left[1+\frac{8}{\lambda ^2}\mathrm{cosh}^2\left(\frac{\lambda }{2}\right)\left(\underset{N=1}{\overset{2n1}{}}{}_{}{}^{}\beta _{N}^{}\right)\right],$$
(37)
where

$$\alpha _N=\frac{1+\frac{4\lambda ^2}{N^2\pi ^2}}{\left(1+\frac{\lambda ^2}{N^2\pi ^2}\right)^2},\beta _N=\frac{5\frac{N^2\pi ^2}{\lambda ^2}}{\left(4+\frac{N^2\pi ^2}{\lambda ^2}\right)\left(1+\frac{N^2\pi ^2}{\lambda ^2}\right)^2},$$
(38)
and we have written the approximation in terms of $`\lambda =gFs`$, $`(\mathrm{\Lambda }_N=\frac{N\pi }{\lambda })`$. That (37) is an approximation to (33) can be seen with the help of the relations

$$\mathrm{cosh}x=\underset{N=1}{\overset{\mathrm{}}{}}{}_{}{}^{}(1+\frac{4x^2}{N^2\pi ^2}),\frac{x}{2}\mathrm{tanh}x=\underset{N=1}{\overset{\mathrm{}}{}}{}_{}{}^{}\frac{1}{1+\frac{N^2\pi ^2}{4x^2}},$$
(39)
whence

$$\underset{N=1}{\overset{\mathrm{}}{}}{}_{}{}^{}\alpha _{N}^{}=\frac{\mathrm{cosh}\lambda }{\mathrm{cosh}^2\left(\frac{\lambda }{2}\right)},\underset{N=1}{\overset{\mathrm{}}{}}{}_{}{}^{}\beta _{N}^{}=\frac{\lambda }{8}\mathrm{tanh}\lambda \frac{\lambda ^2}{8}\text{sech}^2\left(\frac{\lambda }{2}\right).$$
(40)
The nature of the approximation is now evident. The function $`\frac{\mathrm{sinh}\lambda }{\lambda }`$ is rewritten as $`\frac{\mathrm{cosh}\lambda }{\mathrm{cosh}^2\left(\frac{\lambda }{2}\right)}\left\{1+\frac{8}{\lambda ^2}\mathrm{cosh}^2\left(\frac{\lambda }{2}\right)\left[\frac{\lambda }{8}\mathrm{tanh}\lambda \frac{\lambda ^2}{8}\text{sech}^2\left(\frac{\lambda }{2}\right)\right]\right\}`$, the foremost factor is expressed exactly as the infinite product of the $`\alpha _N`$s, and the expression in the square brackets as the infinite sum of the $`\beta _N`$s. The $`n`$th approximation is then defined by including the first $`n`$ terms in the product and in the sum.

Note that $`detA^{(n)}`$ is an even function of $`\lambda `$ (of $`F`$), and hence symmetric, for all $`n`$. This is desirable since the exact (33) is symmetric. The zeroth approximation, in which (33) is approximated by $`detA^{(0)}=2s`$, is the only perturbative result, of order $`(gFs)^0`$. Every approximation order greater than zero contains all (natural number) powers of $`gFs`$.

The matrix (36) shall be our starting point for finding $`detA_1^{(n)}`$. Before replacing the first row of (36) with the column-reduced form of $`C^{(n)T}(s)`$ (the transpose of (35) with $`s^{}=s`$), we must reinstate the second row via row$`{}_{2}{}^{}`$ row$`{}_{2}{}^{}\frac{1}{2}gF`$row1. This procedure is valid since the second element of $`C^{(n)}(s^{})`$ in (22) is proportional to the first element, and we could have used the second element to eliminate the $`e^{gFs^{}}`$ terms instead of the first element. After thus reinstating the second row and replacing the first row with the column-reduced $`C^{(n)T}(s)`$, we obtain the reduced form of the (neither symmetric nor antisymmetric) matrix $`A_1^{(n)}`$:

$$\left[\begin{array}{ccccccccc}e^{gFs}& 0& \frac{1}{2}& 0& \frac{1}{\sqrt{s}}\frac{1}{1+\mathrm{\Lambda }_1^2}& \mathrm{}& 0& \frac{1}{\sqrt{s}}\frac{1}{1+\mathrm{\Lambda }_N^2}& \mathrm{}\\ gFs& 0& \frac{1}{2}(e^{gFs}+1)& \frac{1}{2}gFx_1^T& 0& \mathrm{}& \frac{1}{2}gFx_N^T& 0& \mathrm{}\\ \frac{1}{gF}(e^{gFs}1)& 1& 0& 0& 0& \mathrm{}& 0& 0& \mathrm{}\\ x_1& 0& 0& \mathrm{\Lambda }_1^2& \frac{3\mathrm{\Lambda }_1}{1+\mathrm{\Lambda }_1^2}& \mathrm{}& 0& 0& \mathrm{}\\ 0& 0& 0& \frac{3\mathrm{\Lambda }_1}{1+\mathrm{\Lambda }_1^2}& \frac{5\mathrm{\Lambda }_1^2}{(1+\mathrm{\Lambda }_1^2)^2}& \mathrm{}& 0& 0& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \\ x_N& 0& 0& 0& 0& \mathrm{}& \mathrm{\Lambda }_N^2& \frac{3\mathrm{\Lambda }_N}{1+\mathrm{\Lambda }_N^2}& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}& \frac{3\mathrm{\Lambda }_N}{1+\mathrm{\Lambda }_N^2}& \frac{5\mathrm{\Lambda }_N^2}{(1+\mathrm{\Lambda }_N^2)^2}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right],$$
(41)
We may now expand the determinant to obtain the following $`n`$th approximation to (34):

$$detA_1^{(n)}=\frac{1}{2}\left(\underset{N=1}{\overset{2n1}{}}{}_{}{}^{}\alpha _{N}^{}\right)\left[1+\lambda +e^\lambda +\frac{16}{\lambda }\mathrm{cosh}^2\left(\frac{\lambda }{2}\right)\left(\underset{N=1}{\overset{2n1}{}}{}_{}{}^{}\gamma _{N}^{}\right)\right],$$
(42)
where

$$\gamma _N=\frac{1\frac{2N^2\pi ^2}{\lambda ^2}}{(4+\frac{N^2\pi ^2}{\lambda ^2})(1+\frac{N^2\pi ^2}{\lambda ^2})^2}.$$
(43)
Using the second relation of (39), we find

$$\underset{N=1}{\overset{\mathrm{}}{}}{}_{}{}^{}\gamma _{N}^{}=\frac{\lambda }{8}\mathrm{tanh}\lambda \frac{\lambda }{8}\mathrm{tanh}\left(\frac{\lambda }{2}\right)\frac{\lambda ^2}{16}\text{sech}^2\left(\frac{\lambda }{2}\right).$$
(44)
The function $`2e^\lambda `$ is thus rewritten exactly as

$`{\displaystyle \frac{\mathrm{cosh}\lambda }{\mathrm{cosh}^2\left(\frac{\lambda }{2}\right)}}\left\{1+\lambda +e^\lambda +{\displaystyle \frac{16}{\lambda }}\mathrm{cosh}^2\left({\displaystyle \frac{\lambda }{2}}\right)\left[{\displaystyle \frac{\lambda }{8}}\mathrm{tanh}\lambda {\displaystyle \frac{\lambda }{8}}\mathrm{tanh}\left({\displaystyle \frac{\lambda }{2}}\right){\displaystyle \frac{\lambda ^2}{16}}\text{sech}^2\left({\displaystyle \frac{\lambda }{2}}\right)\right]\right\},`$the first factor again expressed as the infinite product of the $`\alpha _N`$s, the expression in the square brackets as the infinite sum of the $`\gamma _N`$s, and approximated by taking the first $`n`$ terms in the product and in the sum. Note that the zeroth approximation $`detA_1^{(0)}=\frac{1}{2}(1+gFs+e^{gFs})`$ to (34) is non-perturbative.

The $`n`$th non-perturbative approximation to $`\frac{1}{2}(gF\mathrm{coth}gFs+gF)`$ is the ratio of (42) and (37) (note that the product over the $`\alpha _N`$s cancel) is and is the object that appears in the $`n`$th approximation to the Green’s function (32). It is not obvious from the reduced form of the matrix $`A_1^{(n)}`$ that the symmetric part of the determinant of the same is the determinant of the matrix $`\stackrel{~}{A}^{(n)}`$. That is no problem, we find $`det\stackrel{~}{A}^{(n)}`$ in much the same manner as we found $`detA_1^{(n)}`$: take the reduced form of $`A_1^{(n)}`$, reinstate the second column with the operation col$`{}_{2}{}^{}`$ col$`{}_{2}{}^{}+\frac{1}{2}gF`$col1, and discard the first row and column to obtain the symmetric, reduced form of $`\stackrel{~}{A}^{(n)}`$:

$$\stackrel{~}{A}^{(n)}\left[\begin{array}{cccccccc}0& \frac{1}{2}(e^{gFs}+1)& \frac{1}{2}gFx_1^T& 0& \mathrm{}& \frac{1}{2}gFx_N^T& 0& \mathrm{}\\ \frac{1}{2}(e^{gFs}+1)& 0& 0& 0& \mathrm{}& 0& 0& \mathrm{}\\ \frac{1}{2}gFx_1& 0& \mathrm{\Lambda }_1^2& \frac{3\mathrm{\Lambda }_1}{1+\mathrm{\Lambda }_1^2}& \mathrm{}& 0& 0& \mathrm{}\\ 0& 0& \frac{3\mathrm{\Lambda }_1}{1+\mathrm{\Lambda }_1^2}& \frac{5\mathrm{\Lambda }_1^2}{(1+\mathrm{\Lambda }_1^2)^2}& \mathrm{}& 0& 0& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \\ \frac{1}{2}gFx_N& 0& 0& 0& \mathrm{}& \mathrm{\Lambda }_N^2& \frac{3\mathrm{\Lambda }_N}{1+\mathrm{\Lambda }_N^2}& \mathrm{}\\ 0& 0& 0& 0& \mathrm{}& \frac{3\mathrm{\Lambda }_N}{1+\mathrm{\Lambda }_N^2}& \frac{5\mathrm{\Lambda }_N^2}{(1+\mathrm{\Lambda }_N^2)^2}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right],$$
(45)
Expanding the determinant yields

$$det\stackrel{~}{A}^{(n)}=\mathrm{cosh}^2\left(\frac{\lambda }{2}\right)\left(\underset{N=1}{\overset{2n1}{}}{}_{}{}^{}\alpha _{N}^{}\right).$$
(46)
Upon comparison of (46) and (42), noting that $`\alpha _N`$, $`\beta _N`$, and $`\gamma _N`$ are all even functions of $`\lambda `$ ($`F`$), and using $`\frac{1}{2}(1+\mathrm{cosh}\lambda )=\mathrm{cosh}^2\left(\frac{\lambda }{2}\right)`$, our earlier statement that $`det\stackrel{~}{A}^{(n)}`$ is the symmetric part of $`detA_1^{(n)}`$ is apparent.

Finally then, (42), the symmetric part thereof, (46), and (37) substituted into (32) give the $`n`$th approximation to the fermion Green’s function in a constant external field, in the Schwinger-Fock gauge, as defined through the exact representation (4)-(6). The approximation easily leads to Schwinger’s exact result (8) in the limit $`n\mathrm{}`$, and we repeat and emphasise that only the zeroth approximation corresponds to a perturbative result.


text:
# Effective Potential for Polyakov Loops in Lattice QCD

## 1 Introduction

Precise measurement of Polyakov loop correlation is important for finite temperature physics. First of all heavy quark potentials at finite temperature are obtained by measuring two-Polyakov-loop correlation functions. It is already well-known that we can obtain the $`Q\overline{Q}`$ potential below and above the critical temperature ($`T_C`$) from this measurement. Due to the nonperturbative effects, however, there are still less-known quantities at intermediate temperature. For example the $`Q\overline{Q}`$ potential at finite temperature is generally written as $`V/T=e/(RT)^d\mathrm{exp}(\mu R)`$, where $`T`$ is temperature, $`R`$ the distance between quarks, and $`e,d,\mu `$ are parameters. At very high temperature, perturbation theory predicts $`d=2`$, while just above $`T_C`$, lattice simulation gives the deviation from $`d=2`$. In relation to such a nonperturbative behavior at finite temperature, our main interest here is to derive an effective theory for Polyakov loops from lattice QCD. Recently new Polyakov loop model is proposed regarding the deconfining phase transition in $`SU(N)`$ gauge theory. It conjectures a relationship between the Polyakov loops and the pressure. The usual Polyakov loop $`\mathrm{}_1`$ transforms under a global $`Z(N)`$ transformation as $`\mathrm{}_1e^{i\varphi }\mathrm{}_1`$, with $`\varphi =2\pi j/N(j=1\mathrm{}N1)`$. This means that $`\mathrm{}_1`$ transforms as a field with charge one. On the other hand one can consider Polyakov loops with higher charge under a global $`Z(N)`$ transformation. For example the charge two Polyakov loop is defined as

$$\mathrm{}_2=\frac{1}{N}\mathrm{tr}L^2\frac{1}{N^2}(\mathrm{tr}L)^2,$$
(1)
where $`L`$ is an ordered link product in the temporal direction, $`L(x)=\mathrm{P}\mathrm{exp}(ig_0^\beta A_0(x,\tau )𝑑\tau )`$. The Polyakov loops with higher charge might affect non-universal behavior in the model. Therefore it is interesting to measure relevant kind of new Polyakov loop correlation on the lattice without any phenomenological assumptions. In this case, however, we need the Polyakov loops with relatively large length in the temporal direction, which might be inaccessible in the ordinary computation of the Polyakov loops. Thus we are led to try the new algorithm recently developed by Luscher and Weisz called the multi-level method to measure various Polyakov loop correlation functions. This method can measure quite precisely for the Polyakov loop correlation with the large temporal length and/or large distance between the Polyakov loops as shown below.

## 2 Precise measurement of Polyakov loop correlation

One way to get the Polyakov loop correlation functions as precisely as possible is the method of improved estimator. When we measure Polyakov loops, instead of a time-like link variable itself, by using some average of the time-like links, we can reduce the error of the Polyakov loop without affecting the average. The multi-level algorithm, which is used here, is more efficient in measuring the Polyakov loop correlation. We explain this algorithm briefly in the following. Let us consider the two-Polyakov-loop correlation, $`P(n)P^{}(n+r)`$, where $`n`$ is a spatial coordinate and $`r`$ is the distance between two quarks. First we define a two-time-like link operator,

$$T(t_0,r)_{\alpha \beta \gamma \delta }=U_4(t_0,n)_{\alpha \beta }U_4^{}(t_0,n+r)_{\gamma \delta }.$$
(2)
The essence of the multi-level algorithm is similar to the multi-hit one, i.e., is to take an average of $`T`$ instead of $`T`$ itself in the sublattice of the time slice,

$$[T(t_0,r)_{\alpha \beta \gamma \delta }]=\frac{1}{Z_s}D[U]_sT(t_0,r)e^{S[U]_s}.$$
(3)
This is schematically shown in Fig.1. Here $`s`$ denotes the sublattice which includes the operator $`T`$ and is shown as the shaded area in Fig.1.

The Polyakov loop correlation function is then expressed as

$`P(n)P^{}(n+r)`$ (4)
$`=`$ $`[T(0,r)T(a,r)][T(2a,r)T(3a,r)]`$
$`\mathrm{}[T(T2a,r)T(Ta,r)]`$
for single layer or

$`P(n)P^{}(n+r)`$ (5)
$`=`$ $`\left[[T(0,r)T(a,r)][T(2a,r)T(3a,r)]\right]`$
$`\mathrm{}\left[\mathrm{}[T(T2a,r)T(Ta,r)]\right]`$
for two layers. These equations are also schematically shown in Fig.1.

In this method there are essentially two new parameters, which we call $`n_{\mathrm{ms}}`$ and $`n_{\mathrm{up}}`$. $`n_{\mathrm{ms}}`$ is the number of measurement to be made for the time slice average

$`[T(t_0,r)T(t_0+a,r)]`$ (6)
$`=`$ $`{\displaystyle \frac{1}{n_{\mathrm{ms}}}}{\displaystyle \underset{i}{\overset{n_{\mathrm{ms}}}{}}}T_i(t_0,r)T_i(t_0+a,r),`$
and $`n_{\mathrm{up}}`$ is the number of time slice updates between measurement. We determine $`n_{\mathrm{ms}}`$ so that the s/n ratio is equal to unity. Then the statistical errors are reduced exponentially as the distance between the Polyakov loops increases.

We show some test results of the algorithm in the following. All the simulations are done with the lattice size $`12^4`$ at $`\beta =5.7`$. We use the standard Wilson action and no gauge fixing. Fig.2 shows the $`n_{\mathrm{ms}}`$ dependence of a two-Polyakov-loop correlation function at a certain gauge configuration. In the figure, the spatial distance between the two Polyakov loops is set to $`4`$ in the lattice unit and $`n_{\mathrm{up}}=1`$ is fixed. In the single layer case, the relative computation time in the horizontal line is proportional to $`n_{ms}`$, i.e., the leftmost point corresponds to $`n_{\mathrm{ms}}=10`$ and the rightmost $`n_{\mathrm{ms}}=1000`$. In the two layer case, we have chosen the parameters $`n_{\mathrm{ms}}=10`$ and $`n_{\mathrm{up}}=16`$ for the outer layer and $`n_{\mathrm{ms}}=5150`$ and $`n_{\mathrm{up}}=1`$ for the inner layer. Usual measurement of the Polyakov loop correlation without any improvement corresponds to $`n_{\mathrm{ms}}=1`$. The graph shows that the values approach the precise one as $`n_{\mathrm{ms}}`$ increases. We have confirmed that our calculation agrees with that of  within the range of the statistical error.

We have also calculated three-Polyakov-loop correlation shown in Fig.3. As an example of the result we show the $`n_{\mathrm{ms}}`$ dependence of the quantity $`|P(n+4x)P(n+4y)P(n+4z)|`$ in Fig.4. We see the validity of the algorithm even in this case by increasing the value of $`n_{\mathrm{ms}}`$.

Finally we comment on our trial computation of a variant of the multi-level method. Measurement of many-point Polyakov loop correlation in the above method spends much memory. One possible solution is to apply a method which is used in the calculation of quark propagators of disconnected diagrams. Let us consider for simplicity the two-Polyakov-loop correlation. The modified method is to take summation over space of the two-time-link operators before taking the tensor product of them,

$`P(n)P^{}(n+r)={\displaystyle \frac{1}{N_{\mathrm{conf}}}}\{[[K(0)][K(2a)]]`$ (8)
$`\mathrm{}[[K(T4a)][K(T2a)]]\},`$
$`K(t_0)={\displaystyle \underset{\mathrm{space}}{}}T(t_0,r)T(t_0+a,r).`$
If the number of configuration and $`n_{\mathrm{ms}}`$ are sufficiently large, this method gives the same results as the conventional multi-level method, because the gauge non-invariant terms cancel out. Our trial simulation is, however, much noisier than the conventional one due to the finite number of configuration and $`n_{\mathrm{ms}}`$. So unfortunately our present results imply that this method is less effective.

In summary, toward the construction of the effective theory for Polyakov loops, we have tried a new algorithm called the multi-level method. We have reproduced the two-Polyakov-loop correlation by  and confirmed the validity for the computation of the three-Polyakov-loop correlation.