Initial states and apodisation for quantum field simulations in phase-space

Bosonic quantum fields can be simulated with `quantum software' in phase-space. The positive-P, Wigner and Q-function phase-space methods are reviewed. Initial quantum states and boundaries for infinite domains are considered in detail. The quantum initial conditions treated include both multi-mode Gaussian states and number states, which are sampled using phase-shifted weighting techniques. This is more efficient than direct probabilistic sampling. Algorithms for treating periodic boundaries within infinite domains are developed via apodisation, or absorbers near the boundary. Similar truncation or anti-aliasing issues arise in the frequency domain with Fourier transform methods. Complex apodisation, in which there are phase-shifts as well as absorption, improves accuracy by reducing unwanted errors, and this is illustrated numerically. The complex power law apodisation techniques described here can be applied to numerical calculations or even directly in optics. Phase-space simulations, used for both quantum optical and Bose-Einstein condensate quantum dynamical calculations, require quantum noise terms to be included. A method is described to analyze number conservation in apodised calculations. Some techniques developed here are potentially applicable to all digital simulations, whether using classical or quantum logic.


Introduction
Phase-space methods are widely used in quantum simulations. Wigner [1] developed the first quantum phase-space method. Husimi [2] invented a positive representation, and later Moyal [3] derived a dynamical theory. Glauber adapted these approaches to normally-ordered operators in laser physics [4,5]. These techniques were later extended, shown to be equivalent to stochastic processes [6] and used to simulate quantum fields in super-fluorescence, pulsed quantum squeezing [7,8], and many other quantum processes [9]. This paper describes essential components of space-time quantum field dynamical simulations. These have been used to treat very large Hilbert spaces, up to to millions of qubit equivalents [10]. In order to utilize such techniques 1 Introduction one must both generate initial quantum states, and implement boundary conditions. Gaussian initial states are relatively simple to implement. The important case of initial number states is also treated here. For both initial and boundary conditions, techniques using complex weights and phase-shifts are shown to be effective.
Fast Fourier transform (FFT) techniques [11] are frequently used with phasespace representations to numerically propagate quantum fields. Due to their accuracy and speed, FFTs are used in classical optics, laser physics, photonics and BEC calculations. The use of such methods, especially in infinite domains, requires the treatment of truncation in ordinary and momentum space. Other methods of integration, with non-uniform grids and finite differences or related spectral techniques, are also possible [12]. These alternatives do not alleviate the fundamental issue, common to all field simulations, which is the problem of representing a large space-time volume with a finite number of lattice points.
FFT techniques usually employ periodic boundary conditions, which permit an efficient solution in quantum field simulations. Thermal and quantum noise can be easily included. However, periodic boundary conditions can cause errors. This is especially an issue with localized initial fields on infinite backgrounds, where a periodic boundary condition is not physical. The fields can propagate through the boundaries and wrap around. There are similar "wrap-around" or aliasing errors in frequency space. This is a generic problem, and holds not just in the case of transverse diffraction. It also occurs in the case of short pulse propagation, where dispersion effects play a similar role to diffraction. Again, it is usually most efficient to use Fourier transform techniques on a finite window, with periodic boundary conditions. We regard temporal windowing as a space window, just as in the case of transverse diffraction effects.
On long time-scales, periodic boundary conditions in space lead to wraparound that is not present in the target system. Similar issues also arise in optical instruments of many kinds, where a hard boundary causes unwanted diffraction and interference [13]. While such devices do not have periodic boundaries, the use of gradual absorption or apodisation has long been used to alleviate problems due to diffraction at hard boundaries. In the software world, boundary absorption or apodisation is a useful way to prevent boundary artifacts in either momentum or space [14,15].
For quantum software using phase-space methods, when simulating quantum fields, one must also consider the effects of apodisation on noise and conservation laws. Without taking account of this, errors can occur in the field commutation relations and in calculating the total numbers of particles. Any proposal must satisfy the requirement that in a some neighborhood of x = 0, there is a minimal perturbation of the exact quantum solution.
This paper treats both quantum initial conditions and apodisation methods. In both cases, complex weights and phase-shifts are utilized to improve simulation efficiency. Related techniques are also used in lithography and lens design [16]. Here we show how such ideas can be used to simulate quantum field dynamics with a finite digital machine. This issue is likely to occur in all types of quantum simulations, including possible future quantum computers.

Phase-space representations
We start by summarizing the properties of phase-space representations. The positive-P and truncated Wigner-Moyal phase-space methods have been employed to simulate one-dimensional (1D) quantum field dynamics in photonics [17], and in Bose-Einstein condensate (BEC) quantum dynamics in higher dimensions [18]. Both methods were also applied to full simulations of optical fibers, including Raman effects [19].
Later applications included large-scale truncated Wigner simulations of BEC collisions [20], simulated quantum transport in a BEC in an optical lattice [21], detailed comparison of experiment and theory for optical fiber, to below the quantum noise level [22], comparisons of 3D positive-P and Wigner simulations of BEC collisions with 150,000 atoms from first principles [10], proposals for atomic Hong-Ou-Mandel correlations from BEC collisions [23], violations of Bell inequalities [24,25], and many other applications [9].
The starting point is an M −mode Hilbert space with a number-state basis of |n , where n = (n 1, . . . n M ). The Bose operators areâ, whereâ = (â 1, . . .â M ). Extended vectors are defined as − → a = â,â † , so thatâ † j =â j+M . There is a progression of phase-space representations from normally-ordered, to symmetrically ordered, then anti-normally-ordered approaches. The corresponding phase-space distributions have an increasing variance as one goes from normal to anti-normal ordering.

P-representation
The Glauber P-representation [4] uses a diagonal coherent projector, and is singular for non-classical fields. The generalized P-representation [6] extends this by using an off-diagonal coherent projector to treat all quantum states without singularities. It is defined as a real or complex phase-space function P ( α), such that:ρ Here α = (α, β) is a complex vector in a non-classical phase-space of 4M real phase-space dimensions, and dµ (α) is an integration measure. This can either be a volume or surface measure in the enlarged phase-space α. The result is an exact expansion of the quantum density matrix in terms of the coherent state projection operatorΛ This projector is defined in terms of M -mode Bargmann-Glauber coherent states α , where A consequence of the definition is that the P-function is normalized, which follows since T r (ρ) = T r Λ = 1, so that P ( α)dµ (α) = 1. There are existence theorems and explicit constructions for all quantum states, known for both complex and positive forms. The most general case is a complex P ( α), written as where P + ( α) is real, positive and normalized, so it can be sampled with conventional probability techniques: Here Ω( α) is a complex weight, which is ideally chosen so that has a nearly uniform modulus, |Ω( α)| ∼ 1, to give efficient sampling. From these definitions, it has a mean value of unity: It is also convenient to define a weighted mean, which includes the complex weights, so that: Averages in P-distributions correspond quantum-mechanically to normalordered quantum averages, denoted : . . . :. These have creation operators likê a † ordered to the left, and if O( α) is a power series in α, then the following equivalence holds: The average Ô ≡ T r Ôρ is a standard quantum ensemble average with density matrixρ, and O P indicates the P-function average. When sampling the distribution with S samples, averages are obtained from samples of α (n) distributed with probability P + ( α (n) ),so that: The main drawback of this approach is that the sampling error grows in time due to the enlarged phase-space, and in some cases there are boundary term errors in phase-space. This requires the use of additional stabilizing gauges, and limits the time-duration that can be usefully simulated [26,27,28,29].

Wigner representation
Other representations are known, including the Wigner and Q-representations. These correspond to symmetrically ordered and anti-normally ordered operators respectively. The phase-space is classical, with α ≡ [α, α * ]. They can be written formally as convolutions of the generalized P-function, extending a result of Cahill and Glauber [5]: One can define the original Glauber-Sudarshan P-representation as the singular limit of the convolution in the limit of s → 0. However, this often may exist only as a generalized function, depending on the quantum state involved. As a result, we will use the positive-P distribution P ( α) instead, so that all relevant nonsingular cases can be treated.
There are multiple ways to define the Wigner representation. One way to define this is as a convolution of the generalized P-function as given above. Thus, a delta-function P-representation representing a coherent state corresponds to a Gaussian Wigner function.
For the case of the Wigner function one has s = 1/2, so that P 1/2 ( α) ≡ W ( α). The Wigner function is always nonsingular and real. This representation is also generally non-positive, making direct probabilistic sampling impossible except for the special cases described below, or by using weighted sampling methods with non-positive weights.
Averages are defined as integrals over the Wigner function, so if d α is now the complex phase-space volume integral, corresponding to 2M real dimensions, then: Wigner phase-space averages like this correspond quantum-mechanically to symmetricallyordered quantum averages. If we define O(â,â † ) s as an s-ordered average, which could be normally ordered, symmetrically ordered or anti-normally ordered, then: Just as with the P-function, one can also define a weighted version, so that, if W + ( α) is defined as a positive probability, then: Other methods for direct sampling are available in the cases where W (α) is positive, but for pure states these are always Gaussian cases, unless approximations are used. The double-dimension formalism of the positive-P method is also applicable here, if we define a positive Wigner function on an extended phase-space α = [α, β] [30].
Typically, the Wigner dynamical equations are truncated using an 1/N expansion for N bosons, in order to obtain tractable stochastic equations.

Q-function
The Q-function, originally proposed by Husimi [2], is defined either using the convolution method, or directly using the fundamental definition that: As with the Wigner case, the extended vector α is also applicable here if one defines α = [α, α * ]. For the case of the Q-function one has s = 1, so that P 1 ( α) ≡ Q ( α). The Q-function is always nonsingular, real and positive. Direct probabilistic sampling is therefore possible. This method is also applicable to qubits and fermions [? ? ], and has been used for simulating large-order correlations and Bell violations [? ? ].
Averages are defined as integrals over the Q-function, so if d α is the complex phase-space volume integral, corresponding to 2M real dimensions, then: Averages like this correspond quantum-mechanically to anti-normally-ordered quantum averages, ie, While the Q function is always positive, it does not generally have a normal stochastic process. Instead, it corresponds to a stochastic process that propagates simultaneously forward and backwards in time [31,32], requiring more specialized techniques to solve it.

Sampling initial quantum states
To treat large quantum systems, there is a need for computational scalability as the number of modes increases. This means that probabilistic sampling is necessary, as direct high-dimensional integration over a multimode distribution is exponentially complex in the number of modes. Sampling is essential for computability purposes.
Scalability is even an issue when there are exact solutions for the eigenstates, because representing an arbitrary initial quantum state requires a superposition of exponentially many eigenstates, if one has a large number of quantum particles or modes. In practice, this is a severe limitation on dynamical calculations even in exactly soluble models [33].
As a result, computability requires an efficient strategy for sampling a range of suitable initial conditions.

Gaussian states
It is straightforward to sample the phase-space distributions described above if the initial conditions are Gaussian states. Examples of these include coherent states [4], squeezed states [34], Bogoliubov states [21,35,36], and thermal states. In these cases, sampling can be obtained by generating 2M real Gaussian noises defined so that together with corresponding random phase-space samples where α = [α 1 , . . . α 2M ] , such that: Here the quantum mean value is given by: In order to achieve an s-ordered cross-correlation -including the entire extended vector of 2M annihilation and creation operators -of one simply chooses the 2M × 2M complex square root B, such that: There are restrictions on B for the Wigner and Q-functions. The condition that α * j = α M+j implies that: In some "quasi-Gaussian" states, like a laser or a BEC, the true quantum states have no well-defined phase [37]. In such cases, one can generate a Gaussian state relative to some phase-reference amplitude, then randomize the phase. If the dynamics and the final measurements are phase-independent, as is often the case, the last step of phase-averaging is not necessary, as it doesn't change the results. In these cases a coherent state is equivalent to a Poissonian average over number states, and therefore either type of input can be used [36].

Number state sampling
In some cases, one wishes to sample the phase-space for a number state, or a superposition of number states. These are more difficult to sample than Gaussian states, but certainly not impossible. The most general possible quantum states can be characterized by their matrix elements in a number state basis such that: or alternatively, that:ρ = nm C nm |n m| .
For a number state mixture, one would have that all the off-diagonal terms vanish, and:ρ = n C n |n n| .
We treat examples of these in this section.

P-function for number state expansions
While the positive P-distribution P ( α) always exists, this is not always the most numerically efficient choice for sampling. While one can always construct a number state using the positive P-representation [? ], there are other techniques that are more effective in terms of sampling. Choosing the complex P-representation adds a complex weight to the initial condition, while still allowing the usual propagation equations for subsequent time evolution. This method has been used in simulations of boson sampling equivalent to permanents as large as 100 × 100 [38,39,40], with sampling errors that scale better than experiments. It has also been applied to simulations of BEC quantum dynamics in one dimension, with initial quantum number states [41]. This technique was described previously, and we review it in greater detail here.
We focus on the complex P-distribution, which is more compact than the positive P-distribution [6]. The existence theorem states that, for an expansion of an arbitrary density matrix in number states such that C nm = n|ρ |m , and for circular contours enclosing the origin in each amplitude, then by Cauchy's theorem, For vacuum states, the P-function sample is simply one of α k = β k = 0. For non-vacuum inputs we use a circular contour of radius r k and make a transition to polar coordinates, so that α k = r k exp iφ k . The radius r k is chosen to minimize the sampling error, and samples are generated on the contour, to give a Monte-Carlo sampled contour integral.
While the expansion is always applicable, we now assume for simplicity that we can choose the mode basis that defines the mode operators such that the density matrix factorizes withρ = ρ k . For multiple independent modes, a generalized P-representation is defined as: The phase variables can be understood intuitively on defining: where φ k is the classical phase, and θ k is a nonclassical phase which only exists when the quantum state has nonclassical features. Hence, one has the result that: In this case we get the overall density matrix as: where we define C nm (φ k ) = C nm e −i(n k −m k )φ k , andn k = (n k + m k ) /2, so that I nm k = 0 for n k = m k = 0, and otherwise the expansion coefficient I nm k is given by : More general cases with coherences between the modes can be treated, and lead to additional phase-coherence terms in the contour integral sampling.

Contour sampling techniques
To illustrate the general sampling procedure, we suppose that each mode has an input number state, so that if P k is the single-mode distribution for mode k with n k bosons: After separating the real and imaginary parts of the exponential, and transforming to phase variables, one can use random probabilistic sampling for the real part, so that: We call the distribution P k a circular von Mises complex-P distribution (VCP), since the weight on the contour is a von Mises probability distribution [42]. Once the phase angle is randomly chosen, the imaginary part becomes an additional complex weight, Ω k . The sampled probability distribution for the subset of modes with nonzero inputs is then where I 0 (x) is the modified Bessel function of the first kind of order 0, θ k ∈ [−π, π), and φ k ∈ [−π, π). This corresponds to sampling the variables separately as φ k = U (−π, π), and θ k = VM 0, r 2 , where U is the uniform distribution, and VM is the circular von Mises distribution [42]. For each sample we calculate the phase-space variables as α k = r exp (i (φ k + θ k /2)), β k = r exp (i (φ k − θ k /2)). For vacuum modes where n k = 0, we use the deltafunction distribution P k (α k , β k ) = δ (α k ) δ (β k ); that is, one takes α k = β k = 0.
The complex weights associated with each sample are These are used to calculate any moment f (α, β) in conjunction with the S drawn samples of α and β as where Ω = Ω k is the total weight associated with the trajectory.

Large radius limit
It is generally optimal in terms of sampling efficiency to choose r 2 k ≈ n k . For larger particle numbers, these expressions then simplify, because one can use a large radius for the integration contour. The phase distribution becomes asymptotically a standard Gaussian of width 1/r in the limit of large r, and we use the asymptotic Bessel function result that: The sampled probability distribution for the subset of modes with nonzero inputs is then a uniform circular distribution of 1/ (2π), multiplied by a normalized Gaussian: In the same large r limit, the weight can be calculated using Stirling's approximation. Introducing n! = n n √ 2πne −n (1 + O(1/n)) and taking r 2 = n, one then obtains for the k − th complex weighting term: Therefore it follows that Ω = exp −inθ 3 /6 to leading order. Noting that θ 2 = 1/n and hence: nθ 3 ∼ 1/ √ n, these are small phase factors, and the normalization involved is of order: where one can show that θ 6 = 15 n 3 .
As a result, in this approximation the overall weight factor is actually Ω ∼ 1 − 15/ (72n), while it should be Ω = 1. This approximation error of order 1/n can be removed by including higher order terms of O(1/n), if the next order of accuracy is required.

Wigner and Q distribution sampling
There are several ways to sample the Wigner and Q functions, and we will describe two possibilities here. The first method is a direct sampling based on known analytic forms of the function, while the second method is convolution sampling based on the generalized P-distribution as a starting point.
One can obtain an analytic form of the Wigner function in terms of associated Laguerre functions, in the case of a factorizedρ, which is diagonal in the number basis, as: From this result, we can see that while the Wigner function always exists for any quantum state, but it has no generally positive form, even when obtained from a positive P-function. High-order Laguerre functions are difficult to calculate numerically, owing to underflow, overflow and round-off problems. Sampling them is also nontrivial, as they are extremely oscillatory. As a result, this approach is not usually employed for carrying out stochastic integration.
Q-function sampling is more straightforward, as this is always positive. For a number state n, one obtains that: This is essentially a product of gamma distributions, for which efficient sampling methods are known.

Convolution sampling
We will now show that by using a different technique, not based on Laguerre polynomials, that sampling the full Wigner distribution is possible, even for relatively complex quantum states. These techniques have proved useful in establishing the properties of the Wigner function when a full-scale exact simulation of a quantum Schrodinger Cat transfer in optomechanics [43] was carried out using a complex P-distribution, which is an exact method. The results given above for convolutions imply that there is always a general sampling theorem. Both Q( α) and W ( α) can always be obtained using sampling methods provided that an algorithm for sampling the generalized P-function is known, since a convolution is equivalent to adding a random gaussian variable. In the present case one must also take into account that the convolution includes a dimensional reduction from a non-classical phase-space.
To sample either the Wigner or Q-function, we therefore first define the following linear combinations: The variable α + is equivalent to a 'classical' coordinate, while α − indicates the degree of non-classicality. To construct a real Gaussian convolution, if we define The overall distribution is therefore a weighted convolution: where the complex weight function is: Thus, if one samples the P-function to obtain α, then samples the convolution variable, ∆ k = ∆ x k + i∆ y k as real Gaussian variables with variance s/2, such that: the resulting phase-space coordinate of α = α + + ∆ is a Wigner or Q-function sample. However, when calculating ensemble averages with a generalized Pdistribution as the starting point, an additional complex weight must be included with: While this is a complex weight, its average is real.

Space and momentum boundaries
After initializing the quantum state, differential identities are used to obtain dynamical evolution equations in phase-space, which are described in detail elsewhere [44]. These are Fokker-Planck equations in structure, and can either be written in terms of the mode amplitudes α, or the corresponding local fields − → ψ (x). On sampling, the FPEs can be transformed into stochastic partial differential equations (SPDE).
When written in Stratonovich [45] form, the field derivative D − → ψ can be expanded as: Here, A is a vector, B a matrix, L is a linear differential operator in a real space x, and ∆w is a Gaussian noise vector, defined such that: These equations usually involve truncation or other approximations in the Wigner case. The situation is different with the Q-function, which is dynamically equivalent to a forward-backward stochastic equation, requiring more advanced techniques [31,32].
SPDEs are often solved using spectral methods, in which the linear term is converted to an interaction picture [46] using Fourier transforms. Other methods are also possible and are most useful with trapped BECs [12]. Given an , a typical midpoint interaction picture approach using Fourier transforms is to define an operator T = F −1 e L[ik]∆t/2 F ,which propagates the amplitude field with the linear operator by ∆t/2, where F indicates a forward Fourier transform.
The overall algorithm for a step in time starting at − → ψ is therefore a three step process, with a midpoint defined implicitly as − → ψ , and a final result of − → ψ : = T − → ψ (2) .
In calculations of field propagation, apodisation is useful in combination with fast Fourier transform computational methods. It prevents leakage through the spatial boundaries from dispersion or diffraction and similarly prevents leakage in momentum from nonlinear interactions. This would result in unphysical wrap-around of the fields from the positive to the negative coordinate axes, in either space or momentum, due to the use of periodic boundary conditions in space or momentum.
Both in physical optics and in classical field simulations of such equations, apodisation means the addition of an artificial absorber near the boundary. Apodisation in space should be continuous and smoothly varying in order to reduce diffraction artifacts. This is used to smooth or filter the field so that it goes to zero smoothly at the boundary, without additional structure. Because Fourier transforms are symmetrical between position and momentum, apodisation is required both in space and momentum [15].
Just as dispersion leads to a spreading in position space, nonlinearity leads to a spreading in momentum space. There are momentum boundaries at the inverse of the transverse spatial step-size. Such momentum limits cause aliasing errors, where high-momentum components are wrapped around to negative or even low momentum. These delocalized jumps caused by nonlinearities are most effectively de-aliased with projectors that removes high-momentum terms, to give a final linear translation T a = F −1 e L(ik)∆t/2 P (k) F .
A common choice is a projector P (k) such that, for a maximum momentum of k m = π/∆x and a lattice step-size of ∆x: This removes aliasing errors [14,47], which are well-known in many applications of Fourier transforms [48]. More sophisticated algorithms are also possible, that either remove high momentum components more gradually [15], or else require less storage by using a phase-shift technique developed for equations in hydrodynamic studies of turbulence, that cancels only the aliased modes [49].

Complex absorption model
The apodisation behavior required to prevent aliasing in space is different to momentum, because of the different effects of the equations, which produce a localized spread in position from second-order derivative terms. Therefore we consider a generic apodisation as a space dependent absorbing reservoir, which includes a complex absorption coefficient. This modifies the original field equations, so thatψ We will start by considering classical apodisation, and include the effects of quantum noise in subsection (5.4). If applied to a periodic calculation with a transverse lattice spacing of ∆x and boundaries at x m = L/2, we require that there is no absorption at the center, so that γ(0) = 0 , and γ(±x m ) → ∞ at the boundaries. However, it is sufficient to simply have a very high absorption at the boundary, rather than an infinite absorption.
The graphs given below are of the number density ρ i = |ψ i | 2 , where Fig (1) and (2) show the errors occurring in ρ 1 = |ψ 1 | 2 in a non-apodised case. The first figure has a spatial domain of x = ±10. The central error for Gaussian diffraction in a doubled space domain of of x = ±20 is shown in Fig (2). Errors at |x| = 10 are of order of 3% of the initial on-axis intensity despite the enlarged boundaries. This shows that the aliasing errors can only be removed slowly by increasing the spatial volume.
All numerical examples obtained here used a central difference Fourier transform algorithm [46], with public-domain stochastic partial differential equation software [50].
To achieve a minimum error in the central propagation region, the form considered here for apodisation in space is complex, so that it can be expressed in terms of real fields, γ ′ (x) and V (x), so that  Fig. 1: Non-apodised diffraction, Eq (58), in a finite interval with periodic boundary conditions. Graph is of |ψ| 2 versus x, t for a Gaussian of initial width σ = 1, integrated out to t = 20 in a spatial region of ±10.
This complex apodisation, including phase-shifts as well as absorption, is then expanded as a polynomial in the transverse coordinates, where we consider just a single transverse dimension for simplicity: This allows one to investigate a range of possible models of apodisation. Allowing for both real and complex terms corresponds to including arbitrary polynomial potentials V (x) to shape the solution, in addition to the absorption term, γ ′ . These correspond physically to a change in refractive index in optics, or an external potential in quantum mechanics. The challenge is to find an absorption function that minimizes errors caused by the periodic boundaries, and gives a result for r |x| < x m /2 close to the exact one.

Time-evolution
To illustrate this approach, we analyze a relatively simple case in one dimension of the paraxial equation for diffraction or dispersion of a field. To treat this problem, we use an analytic theory. Suppose that there is an evolution equation in an infinite spatial region, which is scaled to have the form of: . Graph is of |ψ| 2 versus x, at t = 20 for a Gaussian of initial width σ = 1, in a spatial region of ±20.
The apodisation method involves adding an absorption γ (t, x), so that: Here the absorption term γ (t, x) is chosen so that the solution to the modified equation can be treated in a finite region with periodic boundary conditions, having a minimum error relative to the exact solution in the central region, in this case for |x| < x m /2.
To analyze the problem, consider a solution in one transverse dimension using a power series expansion for log ψ, where only even terms are assumed. There is an overall negative sign to ensure that coefficients α p are positive for bounded solutions: The expansion coefficients α p (t) must satisfy a recursion relation. As a result, the expansion coefficients evolve in time according to: where the nonlinear term coming from the second derivative is expanded as: Here, β 1 = 2α 2 1 , β 2 = 8α 1 α 2 , β 3 = 12α 1 α 3 + 8α 2 2 , β 4 = (16α 1 α 4 + 24α 2 α 3 ) , and so on. The nonlinear coefficients β q have the universal property that the highest order term appearing for q > 1 is 4qα q α 1 , and all other terms have lower orders in α q .
The solution to the apodisation problem depends on the exact solution that it is intended to approximate. This actually makes it highly nontrivial. Here we apodise an exactly known Gaussian diffraction problem. Such cases have behavior near the boundaries that resembles the unknown solution of interest. Even though the real problem typically has additional nonlinear terms in its timeevolution equation, the extra nonlinear terms are negligible at low intensities near boundaries.
Suppose that only α 0 and α 1 are nonzero without apodisation, so that the exact Gaussian solution has:α This provides a reference solution for the non-apodised case. Initially, let α 1 (0) = 1/ 2σ 2 , where σ is a scale length that defines the size of the Gaussian beam. It is possible to choose an apodisation strategy that is independent of σ, provided there is a large spatial window. With this choice, one obtains: and therefore one obtains the zero-th order solution, as expected, of This is a standard Gaussian diffraction/dispersion solution, which has the property that it spreads out indefinitely. The result is that the central intensity tends to zero at large t, since ψ (t, 0) ∼ 1/ √ t − iσ 2 . If the same equations are solved with periodic boundary conditions, the solution has a conserved total number, and simply becomes a mass of interference fringes as shown in Fig (1). These fringes only decay slowly with increasing boundary size, as shown in Fig  (2).
To treat apodisation, one may consider terms up to p − th order, for some maximum order p, and we define Γ p = γ p x p m as the boundary absorption. We suppose that α 0 , α 1 are nonzero initially. To give a definite example, suppose that α q = 0 initially for q > 1, as in the exact solution given above. The loworder terms are as usual for a Gaussian, and there also must be an envelope α p in the apodised solution, which prevents wrap-around. This is obtained through adding apodizing terms such that α q = 0 for p > q > 1. Fig (3) shows the result of purely absorbing apodisation, with one term of Γ 10 = 10, giving reduced errors of order of 7 × 10 −5 of the initial on-axis intensity.

Complex apodisation
We treat a p-th order apodisation. From matching powers in the two series, there are 1 + p corresponding time-evolution equations.The higher-order terms for q > p are assumed to have relatively small magnitudes, and are neglected. We wish to specify the apodisation so that α 2 and α 0 are unchanged from their free-field behavior given above, to obtain behavior near the center that is invariant under apodisation, with all the other terms zero up to α p . This implies that γ q = 0 for 0 ≤ q < p − 1, and It is essential that α p has a positive real part to apodise the wave-function as it spreads near the boundaries. Hence, γ p−1 should have complex values. This additional quartic potential with V q−1 < 0 provides a linear phase-shift that cancels diffractive shifts induced by the apodisation.
As a result, since the time-evolution of α 1 is known, the p-th term must satisfy:α To obtain an analytic solution, we first consider a constant apodisation γ p , and introduce a new time variable τ = ln t − iσ 2 . The amplitude equation then becomes: On solving this equation, one finds that the p-th order amplitude has a growing behavior in time,which has a well-defined limit after transforming back to the original time variables, such that: Therefore, to prevent possible lower-order amplitude changes, one should include a time-dependent term or order p − 1, such that: At long times, one finds that There is some growth in the truncation error with time, since α p increases linearly. This is not as severe a problem as one might imagine. The growth simply means that there is a truncation near the boundaries, just as one might expect from having an absorbing region. This also leads to an apparent nonconservation of energy. Again, this is physically reasonable, as diffraction leads to escape through the boundaries, which can be monitored numerically as explained below.
In order to understand the scaling with window size, suppose that γ p is defined by the maximum absorption Γ p at the window boundary of x = ±x m , so that The maximum size of the optimal complex term, γ p−1 , reduces as the window size increases, with an asymptotic behavior in time that is purely imaginary and independent of the initial Gaussian width: As a result, a feasible strategy of checking the boundary error is to increase the window size with constant Γ p , where Γ p is the boundary value of the highest order apodisation term. This should result in negligible change in the central   5: Apodised diffraction in a finite interval with periodic boundary conditions, using tenth order apodisation, with Γ 10 = 10, and a complex phase-shift as in Eq (72). Graph is of wave-function errors |∆ψ| = |ψ a − ψ exact | versus x at t = 20. Here the apodised wavefunction is ψ a . All integration parameters as in Fig (2).
region. The required asymptotic form is independent of the scale factor σ, so that this phase-shift technique can be used even when the correct scale factor is only known approximately, as it is for non-Gaussian and nonlinear calculations. Results are shown in Fig(4), with intensity errors reduced to less than 2 × 10 −5 of the initial on-axis intensity, and in Fig (5), which shows that wavefunction errors are reduced to less than 7 × 10 −4 across the entire region of |x| < 10.
A constant apodisation leads to a growing error in α p . It is also possible to consider a time-dependent apodisation. This leads to solutions where α p is asymptotically constant. However, this is not investigated here. Such methods appear to be closely dependent on the exact solution that is targeted. Hence, they may be less useful generally.

Quantum noise terms and number conservation
Because of commutation relations, a quantum apodisation must be accompanied by noise terms, in either the Wigner or Q methods. This is not needed for P-distributions due to normal-ordering. Such terms are also not needed with projective anti-aliasing in momentum space, since the amplitude of the projected high-momentum modes is maintained at zero. In s−ordered simulations, apodizing absorption should be accompanied by noise, so that: where we omit The purpose of the noise is to ensure that the s−ordered quantum fields have a vacuum occupation of s bosons per mode, with Effectively, the apodisation procedure adds an absorbing zero-temperature reservoir that prevents particles from reaching the boundary. Suppose one wishes to record how many particles have been absorbed? In this case, a second, reservoir occupation field is needed, as shown in Fig (6) for the test case of Gaussian diffraction.
The s−ordered quantum apodisation equations (73) have the property that, if we define then the expectation values change in time due to the reservoirs so that: Hence, we can add a reservoir number field that reflects the complementary process in which atoms are transferred to the reservoir while vacuum fluctuations are lost, which satisfies the following Stratonovich stochastic equation: Other parameters as in Fig(4).
As a result, provided there is no other reservoir present, the total particle number is conserved: d dt ρ 1 (x) + ρ 2 (x) = 0.
This allows the total number of particles to be tracked, independent of boundary absorption and apodisation.

Conclusions
This paper discusses "quantum software" applied to the simulation of quantum field propagation with s−ordered phase-space methods. Typically, via probabilistic sampling, one obtains stochastic partial differential equations which can be solved numerically even in Hilbert spaces of large size. There are other, more subtle issues that arise during time-evolution. These are the generation of arbitrary initial quantum states, and the treatment of finite spatial boundaries. For initial Gaussian quantum states, the initial sampling is straightforward. Complex weighting techniques are more efficient when the initial states are number states. This provides a way for sampling distributions for any quantum states with arbitrary ordering. Similarly, we show that complex apodisation methods can provide an efficient means of treating infinite domains using Fourier transform methods. However, in such cases, noise terms are needed to ensure that commutation relations are respected.