Complex Langevin Equations and Schwinger-Dyson Equations
Stationary distributions of complex Langevin equations are shown to be the complexified path integral solutions of the Schwinger-Dyson equations of the associated quantum field theory. Specific examples in zero dimensions and on a lattice are given. Relevance to the study of quantum field theory phase space is discussed.
keywords:Complex Langevin Equation, Schwinger-Dyson Equations
Pacs:11.15.Ha, 02.50.Ey, 05.10.Gg
, Corresponding author.
One frequently comes across complex valued path integral weights in physics. A major example is the study of field theories in Minkowski space. Even in Euclidean field theories the effective action may turn out to be complex, e.g. QCD with a chemical potential. Another case, in which we will be interested in this paper, is the study of the solution space of Schwinger-Dyson equations as the phase space of the associated quantum field theory [1, 2].
Standard Monte Carlo methods do not work in the case, where the path integral weight is not positive-definite. This is the famous “sign problem”. The root of the problem lies in the fact that Monte Carlo methods work on the probabilistic interpretation of the path integral weight , where is the action of the system. In the examples above, the exponent becomes complex and the probabilistic interpretation fails.
Complex Langevin equations were proposed by Parisi  and Klauder  to simulate systems with complex valued path integral weights. For dimensional systems with real weights, a Langevin equation in dimensions may be used to study the partition function of the system. When properly set up, the Langevin process converges to a unique stationary distribution, which is the partition function of the associated system, in the limit of the additional dimension going to infinity. This fact was used by Parisi and Wu in the stochastic quantization of quantum fields . For systems with complex actions, one can still write down a (complex) Langevin equation, as suggested by Parisi and Klauder, but this approach comes with many problems. First of all, it is not certain that the complex Langevin simulation will ever converge to a stationary distribution and if it does, there may be many such stationary distributions, see e.g. [6, 7, 8, 9, 10]. Salcedo noted that these stationary distributions may be constructed by path integrals over contours that connect zeros of the path integral weight . Other authors noted that these stationary distributions satisfy Schwinger-Dyson equations, e.g. [11, 12, 13, 14]. On a completely different track of research, one of the authors with Garcia and Guralnik studied the solution space of Schwinger-Dyson equations and noted that different solutions to Schwinger-Dyson equations may be written as path integrals over contours that connect zeros of the path integral weight [1, 2], exactly as Salcedo suggested for stationary distributions of the complex Langevin equation. Furthermore, the authors of references  and  suggested different solutions of the Schwinger-Dyson equations be interpreted as different phases of the associated quantum field theory and studied this proposal in detail. Salcedo in his paper  made the same suggestion about stationary distributions of the complex Langevin equations in one sentence at the conclusion, but did not elaborate on it. Our aim in this paper is to point out and clarify the connection between these two lines of research and propose @@ complex Langevin equations as a basis of a new numerical method of studying different phases of a quantum field theory. We start by a rather detailed explanation of the mechanism of complex Langevin equations and the problems associated with them in section 2. Section 3 contains the main result of this paper. There, in a zero dimensional setting, we show that the stationary distributions of a complex Langevin equation are the solutions of the Schwinger-Dyson equations for the associated quantum field theory. Furthermore, these solutions may be constructed by changing the integration contour of path integrals from real paths to contours that connect the zeros of the path integral weight on the complex plane. In section 4, we do the trivial generalization of the problem to a lattice and discuss related issues. We conclude by further summarizing the results of  and  and point out the connection between different phases of a quantum field theory and stationary distributions of complex Langevin equations. Based on this observation, we propose complex Langevin equations as a numerical method of studying different phases of a quantum field theory.
2 The Complex Langevin Equation
In this section, we introduce Langevin equations. We discuss zero dimensional quantum field theories with a scalar field for simplicity. Generalization to vector fields and higher dimensions is straightforward, see e.g. . For systems with action , we are interested in calculating expectation values like
using Langevin equations. We will assume to be an analytic function of the complex variable .
Note that our formulation is a Euclidean space formulation, i.e. the path integral weight will be complex only when the action is complex. However in Minkowski space formulations of quantum field theories, the path integral weight will still be complex even though the action is real, since the weight is defined by . When we speak about complex actions in Euclidean space, our results will be applicable to Minkowski space actions after necessary modifications. One has to note that and introduce appropriate factors of ’s in the generating function definition and terms of Schwinger-Dyson equation.
2.1 Real Actions
When the action is real (for the moment is a real field), one can create a stochastic process using a Langevin equation with a unique stationary distribution :
where is a fictitious time and is the Wiener process normalized to satisfy:
Then one can run this Langevin process to calculate the intended expectation values as in equation (1). We first show that this Langevin process really converges to the intended stationary distribution.
Associated with the Langevin process is a probability density
which can be shown (e.g. ) to satisfy the Fokker-Plack equation
Note that the Fokker-Planck equation enables the normalization condition
to be independent of time since
with an appropriate boundary condition on .
Now we look at the asymptotic behavior of as fictitous time goes to infinity. Introducing the quantity
we can rewrite equation (5) as
where the Hermitian operator is called the Fokker-Planck Hamiltonian. Then
is an eigenstate of the Hamiltonian with eigenvalue . Furthermore, (assuming ) it is the ground state since it is nowhere vanishing. Then the time independent eigenvalue equation
has solutions with the property , for . Using these eigenvalues and eigenfunctions, one can write any solution to equation (2.1) as
Since we are looking for solutions that are probability distributions, we apply the normalization condition (6), which sets . This result is obtained by using the orthonormality of the eigenfunctions . Other coefficients are set by the initial probability distribution associated with the random variable at . Then the limit
implies that one recovers the desired stationary probability distribution for the Langevin process as the fictitous time goes to infinity:
Note that this result is independent of the initial conditions. Going back to the expectation value problem (1),
with given by (1) and ergodicity assures the averaging over the path prescription
2.2 Complex Actions
Now we turn to the case where the action is complex. We want to know if we can still use the Langevin equation to calculate desired expectation values. We start by rewriting equation (2) (now being complex) in terms of two real variables and as
where . is again the Wiener process that is normalized to satisfy the mean and variance conditions of equation (3). Note that the equation for has a zero diffusion coefficient (see  for an example where it is not zero), but still is a stochastic equation through its dependence on . Complex Langevin equations may be modified to include a term, the kernel, that may be useful to stabilize the system, e.g [20, 28]. For our purposes, we focus on equation (2.2). Note that we have two different random variables and , therefore the real probability distribution associated with equation (2.2) will be of the form
where we assume to be analytic.
There are two important questions related to this process. The first question is whether the probability distribution converges to a stationary distribution at all,
If it does, does it converge to the desired result:
To understand the difficulties related to the convergence problem, we derive the Fokker-Planck equation. First we note that applying the rules of Ito calculus (see for example ) to the complex Langevin equation (2.2) will give the identity
where the last line follows from the analyticity of . Then using equation (18) one can show that (assuming appropriate differentiability and boundary conditions) satisfies the following Fokker-Planck equation:
A general statement on the existence of a unique zero eigenvalue (stationary) solution for the operator cannot be made. Furthermore, zero eigenvalue solutions may exist in the sense of distributions .
One can assume a complex valued function on the real axis such that
based on the implicit assumption that this equation actually has a solution . The reverse question, existence of a positive given a complex is disscussed in [21, 22, 23]. Using this definition, analyticity of and integration by parts in equation (2.2) one can show that satisfies the pseudo Fokker-Planck equation
which has the same form as that of equation (5).
due to the analyticity of . Inserting this statement into equation (23) and performing the partial integration assuming necessary boundary and differentiability conditions,
Because the action is complex, one cannot make general statements about the spectrum of (see  for the spectral theorem for a limited class of such operators). Furthermore, the relation between the pseudo Fokker-Planck equation and the complex Langevin equation were derived based on certain assumptions. This should be understood in the sense of distributions, meaning only a formal expression of the identity (2.2). Note that is still a stationary solution (i.e. ), but in general the stationary solutions exist as distributions and the uniqueness of stationary solution is not certain .
3 Stationary Distributions of the Complex Langevin Equation and the Boundary Conditions of the Schwinger-Dyson Equation
Despite the difficulties in proving rigorous results, complex Langevin simulations have been used to study many different problems. Interesting cases are those for which the simulation converges to a stationary distribution which is not equivalent to the original complex distribution, e.g. [3, 27, 9, 29, 10, 30]. This must be related to the existence of other stationary distributions. Here we discuss a conjecture related to the stationary distributions by Salcedo  and the well know result that stationary distributions satisfy Schwinger-Dyson identities, e.g. [11, 12, 13, 14]. We show that they follow from one another, based on other work on the boundary conditions of Schwinger-Dyson equations [1, 2]. We will point out the results of these references during our discussion.
We start by assuming that the complex Langevin process has a stationary state. For the stationary distribution, the LHS of equation (2.2) must be zero. Then setting , reproduces the Schwinger-Dyson identities for the Green’s functions of the quantum field theory defined by the action , i.e.
Now define the generating function for the stationary distribution
We will assume that the radius of convergence for this series is nonzero. Then there exists a neighbourhood of around such that the Schwinger-Dyson differential equation holds. i.e.
To see that (29) produces the same identites as (3), substitute the definition of the generating function , equation (28), in the Schwinger-Dyson equation, differentiate with respect to an appropriate number of times and set at the end. One gets the identities (3) order by order at the end of this procedure.
where is a contour over the complex plane. Inserting this into equation (29) one gets:
This equation can be solved for
and is contour that connects the zeros of on the complex plane.
Now consider polynomial actions,
The contours will be defined by wedges, where is the order of , such that as . Contours obtained by deforming without crossing singularities of and keeping boundary points fixed result in the same generating function. Note that set of all will not be independent. The Schwinger-Dyson equation will be of order , and will have independent solutions. For example, if , then will have three zeros in the complex plane, Figure 1,which will define three different generating functions. However the Schwinger-Dyson equation will be a second order linear differential equation, which has two independent solutions. Therefore any two of the three possible paths will define an independent solution set.
Since the Schwinger-Dyson equation is linear, any linear combination of the (independent) solutions will also be a solution,
where define an independent subset of paths . Now defining the distribution on the real plane as
which can always be done by a real parametrization of the complex contour , one sees that the equilibrium distribution can be written as a linear combination
which is exactly the conjecture that was made by Salcedo , where he derived the same result for a general complex distribution by considering the stationary solutions of the pseudo Fokker-Planck equation (24) to be realized as distributions rather than functions. Actually, we managed to refine his conjecture (which considers a sum over all instead of on the RHS of equation (36)) by showing that not all of are independent through the use of Schwinger-Dyson equations. A final note is that the coefficients may depend on initial conditions. We will illustrate these points with numerical examples in the next section.  has other examples discussed along the lines mentioned here.
3.1 Zero Dimensional Examples
We consider two examples here,
In both cases, we will derive an independent set of generating functionals and identify the particular solution of the Schwinger-Dyson equation to which the simulation converges by observing the sampling points in the complex plane. We will see that the change of initial conditions may change the resulting stationary distribution. Zero dimensional field theories have been heavily studied with complex Langevin equations before, e.g.[26, 3, 27, 28, 13, 29, 30, 9]. What makes our presentation different from the previous studies is the relation to complex path integral solutions of Schwinger-Dyson equations.
where is the time step, and is a Gaussian random variable with zero mean and unit variance satisfying
All our simulations run from to with . We start calculating expectation values after . Error bars stand for the standard deviation of runs.
For , an independent set of solutions to the Schwinger-Dyson equation (29) can be written by connecting the three zeros shown in Figure 1, i.e. , and . We choose the following generating functions and :
We expect the result of the complex Langevin simulation to converge to a linear combination of the distributions defined by these generating functions.
Table 1 shows the results of simulations for this theory. We compare with the exact results for the Green’s functions of and . Figure 2 shows a sample path for the simulation of this action. We repeated the simulations with different initial conditions, some of which are given in Table 1. For the converging simulations, the sample paths localized around the same region as of Figure 2 and the obtained numerical values were similar to those of Table 1. For some initial conditions with large values (e.g. ), we observed nonconverging paths, which could be made to converge (and localized in the same region of Figure 2) by decreasing the step size. This behavior can be understood by inspecting the deterministic part of the complex Langevin equation (i.e. without a noise term in equation (2.2)). The solution to to the deterministic part will be:
where is the initial condition. We see that is a global attractor for all points except the positive imaginary axis. Any path starting from the positive imaginary axis will go to infinity staying on the imaginary axis (i.e. ) in finite time. When the noise term is included, which points along the real axis, these diverging paths will come out the positive imaginary axis and eventually approach the sampling region shown in Figure 2. However, numerically these paths may cause a problem. When the step size is not small enough, the simulation may go to infinity around these points in finite time. This is called the “runaway solution” problem, see  for other examples and more details. Despite this numerical problem, which can be cured by smaller step sizes, the simulations suggest that for the action the complex Langevin algorithm always converges. Inspecting Table 1 we see that the distribution defined by has correlators within the error range of numerical data. The sample path of Figure 2 shows that the simulation does sample around the path of . It definitely does not sample around the positive imaginary axis. Based on these observations, we conjecture that the complex Langevin simulation for theory always converge to the distribution defined by the generating function .
|no convergence||no convergence||no convergence||no convergence|
Now we turn to . Note that this action is real, however does not define a probability distribution on the real line, as the integral is divergent. The complex evolution of the associated Langevin equation is introduced by a choice of complex initial conditions. This procedure, choosing complex initial conditions, will turn every real Langevin process to a complex Langevin process. The statements we made about complex Langevin equations are valid for these cases also.
An independent set of generating functions for this theory can be written by connecting the four zeros of on the complex plane, i.e. , , and . We choose the following generating functions , and :
We expect to see the complex Langevin simulation converge to a linear combination of the distributions defined by these generating functions.
Table 2 shows the simulation results with different initial conditions compared with the correlators of the generating functions. In contrast to the previous case, we see that the results of the simulation is initial value dependent. Initial points on the real axis do not converge at all. Initial values above and below the real axis converge, but the sample points localize in different regions of the complex plane giving different results, see Figures 3 and 4. We can understand this behavior again by inspecting the deterministic part of the complex Langevin equation. This time the solution will be:
where is the initial condition and the square root function gives the principal root. We see that is a global attractor for every point on the complex plane except the real line. The origin will repell any path starting on the real line; these paths will diverge in finite time. When the real noise term is added, simulations with real initial points will stay on the real line and due to the repulsion they will diverge. Initial points on the upper/lower half of the complex plane will be attracted by the origin and the simulations will sample in the upper/lower half plane as seen in Figures 3 and 4. As a result, upper half plane initial points will converge to a different solution of the Schwinger-Dyson equation than those of the lower half plane initial points. Table 2 shows two simulations starting with and . We see that the correlators of are in the error range of the former initial condition, while the correlators of are in the range of the latter simulation. Based on these observations, we conjecture that for , the complex Langevin equation will diverge if the initial condition is real, converge to the distribution defined by if the initial condition is on the upper half of the complex plane or converge to the distribution defined by otherwise.
One might suspect that either or both of the two stationary states we discussed are quasi-stationary states and consequently expect to see the system converge to a unique stable stationary state after a long enough simulation time. Assuming convergence to a stationary state, we will argue that this is not the case, initial points above and below the real line will behave differently in the whole range of simulation. Since the path of the complex field is continuous (but not differentiable), a complex Langevin process starting from above the real line will never end up below the real line. The reverse statement holds also. The reason is that any path going from one half plane to the other must pass through the real line, and once the path is on the real line, it stays on the real line. Both the noise term and the drift term points along the real line. Furthermore, the path on the real line will show nonconvergent behavior as discussed above. So either, all paths diverge at the end, or, assuming convergence, initial point above and below the real line end up sampling in different regions of the complex plane, converging to different stationary distributions. This initial value dependence is crucial in complex Langevin equations.
Actions and were studied in the context of -symmetric quantum field theories in , where only one of the path integral solutions to the Schwinger-Dyson equation is considered, e.g. . In fact, the authors of , Bernard and Savage, provided a formal proof that the complex Langevin equation for should always converge to the distribution defined by , regardless of the initial condition, which we have also observed in our simulations. The formal proof is based on the methods of  and involves defining of equation (23) as a projection of to the complex integration contour (after necessary partial integrations), as opposed to the real line as of equation (26). Then the associated Fokker-Planck Hamiltonian will have a real and positive spectrum with one zero eigenvalue. Bernard and Savage argued that with slight modifications the same reasoning applies for with the integration contour of . They also noted that in simulations of one has to choose the initial point to be in the lower half of the complex plane or else the numerical simulations will be unstable. In our studies, we observed instabilities for initial points on the real line, which we interpreted to be nonconvergent behavior. Initial points in the lower and upper halves of the complex plane were observed to converge to different probability distributions (solutions of the Schwinger-Dyson equation). Our results suggest that the proof given by Bernard and Savage should be modified to include the effects of different initial conditions and all possible solutions of the Schwinger-Dyson equation. We believe that more attention must be paid to the boundary conditions of when doing the projection.
Another point to note is that in both cases, we could not recover the whole solution set of the Schwinger-Dyson equations from the complex Langevin equation. In general, theories will have independent solutions to the Schwinger-Dyson equations. In the cases that we studied, we could recover only of them. We may need to consider other stochastic processes to recover the whole solution set.
4 Lattice Study
It is trivial to generalize this discussion to a lattice. In particular, consider a general Euclidean scalar field theory with a polynomial potential term,
For simplicity we consider a lattice in one dimension, which we call time and denote by . The generating functional of the theory is now a function of variables, where is the number of lattice points. We denote the generating functional by . is the source at lattice point. We will use to mean . By definition, correlation functions are given by
Next, we introduce the difference operators on the lattice,
where stands for the field value at lattice point. The Schwinger-Dyson equation for the scalar field theory (44) can be written by a coupled set of partial differential equations. Using a centered discretization for time derivatives, the system is composed of a partial differential equation for every lattice point, given by
where the lattice spacing is set to one.
This discretization requires one to set boundary conditions on at initial and final lattice points. There are many possible choices, we choose periodic boundary conditions. We assume that the whole space is filled with a lattice of period .
In the absence of the kinetic term, each lattice point acts independently. Therefore, study of the complex Langevin equation for zero dimensional case immediately tells us what will happen on the lattice. When the kinetic term is included, couplings between lattice points take action and problem is more complicated.
A general solution to the lattice Schwinger-Dyson equation can be written as
where denote contours that connect the zeros of the integrand on complex plane. We again normalize so that . For each lattice point there are independent solutions, which leads to independent solutions to the whole lattice Schwinger-Dyson equation. Because of linearity any combination of these solutions is also a solution. We note that there is no ambiguity in the solution written in this form, one can do the integrations in any order.
We turn to the complex Langevin equation for this problem. We again introduce a fictitious time coordinate . The complex Langevin system is now written in terms of stochastic variables . For each lattice point , there is a stochastic equation:
where are independent Wiener processes normalized as before. A repetition of the analysis of the previous section is sufficient to conclude that the stationary distributions of this set of stochastic equations will satisfy lattice Schwinger-Dyson equations (47).
Let’s consider again the theory, this time in one dimension. The Lagrangian is given by
We again note that this theory is bottomless, a normal (real line contour) path integral solution to the Schwinger-Dyson equations does not exist. However, complex contour contour solutions do exist. Table 3 shows the results of numerical simulations for this theory on a one dimensional lattice. As in zero dimensional case, we see that the complex Langevin equation has at least two different stationary distributions. Choice of initial conditions can alter the stationary distributions. We note that this theory was also studied in  and there initial conditions were restricted to lower half of the complex plane. This led to the observation of only one of the stationary distributions, which was concluded to be the stationary distribution that led to a -symmetric theory. This could be described by choosing of the previous section at each lattice point as the integration contour for the generating functional (4). Table 3 suggests that the other solution is given by choosing at every point.
More studies on the lattice must be done to understand the convergence behavior of complex Langevin equations. Here, as well as initial conditions, boundary conditions on the lattice may also take effect. Some other specific questions are listed in the next section. We leave detailed lattice studies to a coming paper. The point of this section is to demonstrate that the results of the previous section applies to the lattice as well.
5 Discussion and Conclusion
Salcedo  suggested that stationary distributions of the complex Langevin equation may be interpreted as different phases of the associated quantum field theory. Some authors used Langevin and complex Langevin equation to give physical meaning to bottomless actions, e.g. [33, 28, 34, 35, 13]111Some of these references use nonconstant kernels in the complex Langevin equation which enlarges the set of stationary distributions  and changes the Schwinger-Dyson equation, e.g. . Our discussion shows that a rephrasal of these questions is to ask which solutions of the Schwinger-Dyson equations define a phase of quantum field theory. Here we discuss this point, and the relation of different phases to the boundary conditions of Schwinger-Dyson equations [1, 2].
Schwinger-Dyson equations are differential equations and admit more than one solution. Therefore it is necessary to set boundary conditions to specify the particular solution one is looking for. If one is solving for a quantum field theory using Schwinger-Dyson equations, it seems reasonable to choose the boundary condition so that the solution is the standard path integral over real fields. However, in many cases this solution actually will not be the physical one, e.g. symmetry breaking phases. Also, in theories with actions unbounded below, the integrals over real fields are not even convergent. In these cases, it is reasonable to look at other solutions of the Schwinger-Dyson equations and study them as possible generating functionals of the associated quantum field theory. Of course, different solutions require specification of different boundary conditions. One way of specifying different solutions is to consider path integrals over complex paths (as opposed to real paths) that connect zeros of the partition function on the complex plane, as was demonstrated in equation (31). Some of these also happen to be the stationary distributions of the associated complex Langevin equation constructed by Salcedo  as shown in this paper. The problem with this approach is the large number of different boundary conditions/solutions and if all these different solutions define a phase or vacuum of the associated quantum field theory. A possible reduction of the solution set comes from taking the thermodynamic limit of the lattice. These issues are discussed in detail in [1, 2].
In the light of the discussion above, one concludes that to study the phase structure of quantum field theories, one needs to study the different solutions to Schwinger-Dyson equations. This task requires new numerical methods in the study of quantum field theory. One suggestion is the Source-Galerkin method, see e.g.  and references therein. This method proposes an expansion of the generating functional in polynomials of the source term and optimizes this expansion by a Galerkin procedure using the Schwinger-Dyson equation. It is successful in many problems, but also proved to be very difficult in many other cases. Another approach is given by mollification of the path integral weight . The connection between complex Langevin and Schwinger-Dyson equations suggests the use of complex Langevin simulations as a numerical method to study the phase structure of quantum field theories.
Some questions and speculations in this quest are:
Is it possible to know a priori if the complex Langevin equation will converge?
What is the exact relation between the initial condition and the stationary distributions of the complex Langevin equation?
Is it possible to recover the whole solution set of the Schwinger-Dyson equation using the complex Langevin equation? If not, can we use other stochastic systems to recover the whole set? Also, what is special about the recovered solutions?
How are these results modified in the continuum? There are an infinite number of solutions to Schwinger-Dyson equations on the lattice. The continuum limit may cause a collapse in the solution set. This is definitely the case if space time translational invariance is required. [1, 2]. Can one see this collapse using complex Langevin equations? This will be very important in understanding the phase structure of quantum field theories.
We address some of these problems in future work.
The authors would like to thank D. D. Ferrante and D. Obeid for useful discussions and conversations. G. Guralnik would like to thank Z. Guralnik for many conversations related to this work. This work is supported in part by funds provided by the US Department of Energy (DoE) under DE-FG02-91ER40688-TaskD.
-  S. Garcia, G. Guralnik, Z. Guralnik, Theta Vacua and Boundary Conditions of the Schwinger Dyson Equations, hep-th/9612079.
-  G. Guralnik, Z. Guralnik, Complexified path integrals and the phases of quantum field theory, arXiv:0710.1256.
-  G. Parisi, On complex probabilities, Phys. Lett. B 131 (1983) 393–395.
-  J. R. Klauder, Stochastic quantization, in: H. Mitter, C. B. Lang (Eds.), Recent Developments in High-Energy Physics, Springer-Verlag, Wien, 1983, p. 251.
-  G. Parisi, Y. Wu, Perturbation theory without gauge fixing, Sci. Sin. 24 (4) (1981) 483.
-  H. Gausterer, S. Lee, The Mechanism of Complex Langevin Simulations, J. Stat. Phys. 73 (1993) 147–157.
-  H. Gausterer, On the correct convergence of complex Langevin simulations for polynomial actions , J. Phys. A 27 (1994) 1325–1330.
-  H. Gausterer, Complex Langevin: a numerical method?, Nuc. Phys. A 642 (1998) 239–250.
-  L. L. Salcedo, Spurious solutions of the complex langevin equation, Phys. Lett. B 305 (1993) 125–130.
-  S. Lee, The convergence of complex Langevin simulations, Nuc. Phys. B 413 (1994) 827–848.
-  J. Berges, S. Borsanyi, D. Sexty, I. O. Stamatescu, Lattice simulations of real-time quantum fields, Phys. Rev. D 75 (4) (2007) 045007.
-  G. Aldazabal, A. Gonzalez-Arroyo, N. Parga, The stochastic quantisation of u(n) and su(n) lattice gauge theory and langevin equations for the wilson loops, J. Phys. A: Math. Gen. 18 (15) (1985) 2975–2993.
-  M. Ito, K. Morita, Note on Stochastic Quantization of Field Theories with Bottomless Actions, Prog. Theor. Phys. 90 (1993) 187–200.
-  S. S. Xue, The fokker-planck equations in lattice gauge theories, Phys. Lett. B 180 (1986) 275–280.
-  M. Namiki, Stochastic Quantization, Springer-Verlag, Berlin Heidelberg, 1992.
-  C. W. Gardiner, Handbook of Stochastic Methods, Third Edition, Springer-Verlag, Berlin Heidelberg, 2004.
-  P. H. Damgaard, H. Huffel, Stochastic quantization, Phys. Rep. 152 (5&6) (1987) 227.
-  M. Namiki, Basic Ideas of Stochastic Quantization, Prog. Theor. Phys. Suppl. 111 (1993) 1–41.
-  H. Nakazato, Y. Yamanaka, Minkowski stochastic quantization, Phys. Rev. D 34 (2) (1986) 492–496.
-  B. Söderberg, On the complex langevin equation, Nuc. Phys. B 295 (1988) 396–408.
-  L. L. Salcedo, Representation of complex probabilities, J. Math. Phys. 38 (1997) 1710–1722.
-  D. Weingarten, Complex Probabilities on as Real Probabilities on and an Application to Path Integrals, Phys. Rev. Lett. 89 (24) (2002) 240201–+.
-  L. L. Salcedo, The existence of positive representations for complex weights, J. Phys. A: Math. Gen. 40 (2007) 9399–9412.
-  H. Nakazato, Thermal equilibrium minkowski stochastic quantization, Prog. Theor. Phys. 77 (1) (1987) 20–25.
-  J. R. Klauder, W. P. Petersen, Spectrum of certain non-self-adjoint operators and solutions of langevin equations with complex drift, J. Stat. Phys 39 (1985) 53–72.
-  C. Bernard, V. M. Savage, Numerical simulations of pt-symmetric quantum field theories, Phys. Rev. D 64 (8) (2001) 085010.
-  H. Q. Lin, J. E. Hirsch, Monte carlo versus langevin methods for nonpositive definite weights, Phys. Rev. B 34 (3) (1986) 1964–1967.
-  H. Okamoto, K. Okano, L. Schülke, S. Tanaka, The role of a kernel in complex langevin systems, Nuc. Phys. B 324 (1989) 684–714.
-  H. W. Hamber, H. C. Ren, Complex probabilities and the Langevin equation, Phys. Lett. B 159 (1985) 330–334.
-  J. Flower, S. W. Otto, S. Callahan, Complex langevin equations and lattice gauge theory, Phys. Rev. D 34 (2) (1986) 598–604.
-  J. Ambjørn, S. K. Yang, Numerical problems in applying the langevin equation to complex effective actions, Phys. Lett. B 165 (1985) 140–146.
-  C. M. Bender, D. C. Brody, J. H. Chen, H. F. Jones, K. A. Milton, M. C. Ogilvie, Equivalence of a complex PT-symmetric quartic Hamiltonian and a Hermitian quartic Hamiltonian with an anomaly, Phys. Rev. D 74 (2) (2006) 025016.
-  J. Greensite, M. B. Halpern, Stabilizing bottomless action theories, Nuc. Phys. B 242 (1984) 167–188.
-  S. Tanaka, M. Namiki, I. Ohba, M. Mizutani, N. Komoike, M. Kanenaga, Stochastic quantization of bottomless systems based on a kerneled Langevin equation, Phys. Let. B 288 (1992) 129–139.
-  S. Tanaka, I. Ohba, M. Namiki, M. Mizutani, N. Komoike, M. Kanenaga, Stabilization of -Model Based on Stochastic Quantization Method with Kerneled Langevin Equation, Prog. Theor. Phys. 89 (1993) 187–196.
-  D. Petrov, Novel approaches to numerical solutions to quantum field theories, Ph.D. thesis, Brown University (2005).
-  D. D. Ferrante, G. S. Guralnik, Mollifying quantum field theory or lattice qft in minkowski spacetime and symmetry breaking, arXiv:hep-lat/0602013.