The sign problem in density matrix quantum Monte Carlo

Density matrix quantum Monte Carlo (DMQMC) is a recently-developed method for stochastically sampling the $N$-particle thermal density matrix to obtain exact-on-average energies for model and \emph{ab initio} systems. We report a systematic numerical study of the sign problem in DMQMC based on simulations of atomic and molecular systems. In DMQMC, the density matrix is written in an outer product basis of Slater determinants and has a size of space which is the square of the number of Slater determinants. In principle this means DMQMC needs to sample a space which scales in the system size, $N$, as $\mathcal{O}[(\exp(N))^2]$. In practice, there is a system-dependent critical walker population ($N_c$) which must be exceeded in order to remove the sign problem, and this imposes limitations by way of storage and computer time. We establish that $N_c$ for DMQMC is the square of $N_c$ for FCIQMC. By contrast, the minimum $N_c$ in the interaction picture modification of DMQMC (IP-DMQMC) only is directly proportionate to the $N_c$ for FCIQMC. We find that this comes from the asymmetric propagation of IP-DMQMC compared to the symmetric propagation of canonical DMQMC. An asymmetric mode of propagation is prohibitively expensive for DMQMC because it has a much greater stochastic error. Finally, we find that the equivalence between IP-DMQMC and FCIQMC seems to extend to the initiator approximation, which is often required to study larger basis sets and other systems. This suggests IP-DMQMC offers a way to ameliorate the cost of moving between a Slater determinant space and an outer product basis.


I. INTRODUCTION
In a recent study, we showed that the density matrix quantum Monte Carlo (DMQMC) method could be applied to molecular systems, extending it beyond original applications to model systems in condensed matter physics. 1 The use of finite temperature electronic structure methods are becoming increasingly important in applications such as plasmonic catalysis, 2,3 the study of planetary interiors, 4 and solid-state materials 5 , where the temperature dependence is key in obtaining physical and chemical properties, such as phase diagrams and excitation energies. The inclusion of temperature in quantum chemistry methods is difficult because at finite temperatures, more than one state is often occupied, increasing the difficulty of solving the Schrodinger equation. DMQMC joins a growing set of methods including other quantum Monte Carlo methods, 6-11 many body theories [12][13][14] and others in attempting to solve the finite temperature problem that has attracted recent attention amongst quantum chemists. [15][16][17][18][19] Many of these methods, like DMQMC, continue to undergo development. [20][21][22][23][24][25][26][27][28] Widespread adoption of all methods in the FCIQMC family, including DMQMC, is hindered, in part, due to the sign problem. In FCIQMC-based methods, coefficients in the wavefunction (or density matrix, in DMQMC) are sampled by a distribution of walkers. The original FCIQMC paper found that simulations that exceeded a critical walker population were able to successfully resolve the signs of the wavefunction and generate an energy estimate that was exact-on-average; it was a) These authors contributed equally to this paper b) Electronic mail: james-shepherd@uiowa.edu not possible to find accurate estimates from populations lower than the plateau. 29 In FCIQMC, walkers arriving at the same site can be exactly annihilated due to having the discrete basis set; this contrasts a continuous realspace basis where the same approach can be much more difficult. While the wavefunction is still being sampled exactly on average, the signal-to-noise ratio is extremely low and prevents exact estimates from being extracted. A simulation with a growing walker population will have its growth briefly stall out, forming a plateau in the total walker population (N w ) as a function of the simulation iteration, known as the "annihilation plateau", as the simulation establishes the sign of critical elements of the wavefunction. When the population has grown above the plateau, the sign problem is resolved and exact energies can be straight-forwardly collected.
The sign problem in FCIQMC was discussed in depth in the early developmental papers 29,30 for the method before subsequently being systematically studied by Spencer et al., 31 whose work we refer to throughout. This work established the origin of the sign problem as an unphysical Hamiltonian (whose solution does not have a sign problem) and that is unavoidably encountered in undersampled dynamics. There are also attempts to leverage this understanding directly using a fixed-node or trial wavefunction approach. 28,32 The development of the initiator approach in FCIQMC removed the annihilation plateau at a cost of introducing a small error in the energy (removed by increasing the number of walkers in the simulation). 33 The motivation for and derivation of this approximation was related to the alleviation of the sign problem and allowed for a much broader scope of applications. The development of the initiator approximation in DMQMC achieved a similar outcome allowing for our previous work on the uniform electron gas and ab ini-tio molecular systems. 1,34 Subsequently, there were also a wide variety of FCIQMC or FCIQMC-like methods development which are beyond the scope of this work to review in detail. [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50] Large scale implementations of the FCIQMC method and related methods have also been developed and these papers review current challenges and developments for the interested reader. 51,52 Here, we conduct a systematic investigation of the sign problem in density matrix quantum Monte Carlo (DMQMC). We find that the annihilation plateau comes from the same unphysical Hamiltonian as in FCIQMC. We measure these critical walker populations for a test set from the FCIQMC literature and find that DMQMC plateau heights are proportional to the square of the FCIQMC plateau height. However, we also show that by moving to the interaction picture (IP-DMQMC) the plateau heights scale linearly with the FCIQMC plateau heights. Despite being able to control the sign problem, IP-DMQMC can have an issue with the trace population as there is no global estimator for the energy in DMQMC unlike in FCIQMC. To address the collapse of the trace population that occurs even when the sign problem is overcome, we examine the initiator adaptation, showing that it has similar performance as the initiator adaptation in ground-state calculations using FCIQMC.
We find that the reason that IP-DMQMC has this plateau height reduction is that the propagation is asymmetric. Comparing asymmetric DMQMC to IP-DMQMC, we find that the critical populations are the same when a shift is used in DMQMC. While asymmetric DMQMC does appear to have cost savings in the required population compared to symmetric DMQMC, this is offset by the need to sample over the rows (or columns), or, equivalently, β-loops. We believe this shows IP-DMQMC is as effective at solving for finite-temperature energies as FCIQMC is at solving zero-temperature energies. We see this work as complementary to our previous and future studies which develop and apply DMQMC as well as the related work of Rubenstein et al. discussing the sign problem for finite-temperature auxiliary field quantum Monte Carlo. 6,7,53

II. METHODS
In this section, we provide a summary of the methods used here. We begin with the three methods primarily used in this work: DMQMC, interaction picture DMQMC and FCIQMC. We then describe the initiator adaptation. We note now that Hartree atomic units are used throughout.

A. Density matrix quantum Monte Carlo
We begin with the original formulation of DMQMC. 54 Starting with the unnormalized thermal density matrix ρ = e −βĤ (1) whereĤ is the Hamiltonian operator and β = (k B T ) −1 , we can show that the density matrix satisfies the symmeterized Bloch equation by differentiatingρ(β) with respect to β. A Euler update scheme, or finite difference approach, with a finite time step, ∆β, can then be used to find the density matrix at any β, followinĝ We then rewrite Eq. (3) in a basis of outer products of Slater determinants to obtain a matrix form that can be solved stochastically, by evolving a population of particles through the inverse temperature regime. The result is where T ij = −(H ij − Sδ ij ) is the update matrix, and S is a variable shift for population control of the particles in the simulation, explained later in this section. The matrix elements ρ ij = D i |ρ|D j are represented by particles in the simulation, where |D i are Slater determinants in the defined finite basis set. The i and j indices begin at i = 0 and j = 0. A population of particles, N w is then used to sample elements of the density matrix by evolving with respect to β, according to Eq. (4). Integer weights were used in the original FCIQMC algorithm (which DMQMC is based off of) and because we want to make comparison with previous results, this is what we use here.
The simulation starts at β = 0, where the density matrix is the identity matrix. The simulation is then propagated to the desired value of β. At each step, the population is updated following rules for spawning and death of particles, summarized below, while particles of opposite signs on each matrix element are annihilated; these steps are closely analogous to FCIQMC. 29 There are three rules for evolving particles that can described as follows: • Spawning: occurs from one matrix element (ρ ik ) to another (ρ ij ), along both the rows and columns.
• Cloning and death: occur on single matrix elements only, and are designed to increase and decrease the population respectively.
• Annihilation: particles of opposite signs on single matrix elements are removed from the simulation.
Spawning will occur with the probability p s (ik → ij) = ∆β|T kj | 2 and the sign will correspond to sign(ρ ij ) = sign(ρ ik )×sign(T kj ). The same equations will hold for spawning from ρ kj to ρ ij .
The sign of the newly spawned particle is important because as can be seen, the sign of the new particle depends on both the sign of the matrix element where the particle spawned from and the update matrix, T , connecting the two elements. Because of this, the signs of newly spawned particles will not always be sign-coherent, resulting in the manifestation of the sign problem. In order to resolve the signs of the particles on the matrix elements, a system dependent number of particles is required. This will be explored further throughout this work.
Cloning and death occur with a probability given by The population increases if sign(T ii + T jj )×sign(ρ ij ) > 0, and decreases otherwise. Annihilation also occurs on single matrix elements and is used to control the sign problem and particle growth within the simulation, and has been show to be key in overcoming the sign problem. 29,31 A population control must be used, so we introduce a variable shift parameter, that is controlled by The shift update is dependent on N w (β), the total number of walkers at β, A, the number of ∆β steps between shift updates, and ξ, a shift damping parameter. The steps outlined above are repeated until the desired inverse temperature is reached. To obtain estimates of thermodynamic quantities, one averages over many independent simulations, termed "β loops". Then, to find the energies the following expression is used, Ĥ = T r(ρĤ)/T r(ρ), where the numerator and denominator of this equation are sampled separately over the course of the propagation through β, and averaged over the desired amount of β loops. In this work, we solely use the projected estimator and not the shift estimator (because the shift estimator does not converge to the finite temperature energy in DMQMC 54 ).

B. Interaction Picture DMQMC (IP-DMQMC)
The interaction picture variant of DMQMC (IP-DMQMC throughout) was developed to overcome two sampling issues present in the original DMQMC method: the initial density matrix rarely contains the important determinants and the distribution of weight fluctuates rapidly as a function of β. 55 Replacing the density matrix with an auxiliary matrix,f , means that the simulation can be started at a non-interacting density matrix, e −βĤ 0 , rather than the identity, providing a good first approximation to the fully interacting density matrix for weakly-correlated systems. The auxiliary matrix can be written:f whereĤ =Ĥ 0 +V andĤ 0 is a mean-field Hamiltonian.
In this work, we use the Hartree-Fock Hamiltonian for H 0 , though it is possible to use a more general mean-field Hamiltonian. In practice,Ĥ 0 only has diagonal matrix elements in a Slater determinant basis, and e −(β−τ )Ĥ 0 only has diagonal matrix elements at any temperature. It is important to note that this matrix evolves from e −βĤ 0 at τ = 0 to e −βĤ =ρ(β) at τ = β, which means IP-DMQMC only samples the correct distribution at τ = β, so separate simulations are required for each β value. We can differentiate the matrixf with respect to τ to find: This equation can be simulated using the rules above, with one change: the cloning/death probability in the second rule changes to p d (ij) = ∆τ |H 0 ii − H jj |, becausê H 0 is diagonal in the chosen basis. The condition on increasing or decreasing the population is then based on the sign of (H 0 ii −H jj )×sign(ρ ij ), with population increasing if this expression is greater than 0. In this work, IP-DMQMC uses the asymmetric spawning mode as described in Sec. II C below, meaning that spawning is restricted to occur only along rows.
IP-DMQMC is the same as DMQMC, in that many simulations need to be averaged to obtain estimates for observables. When introduced, it was said that one major benefit of this variant is that as long as H 0 ii > H jj , there is little to no death along the diagonal; this overcomes one problem with large systems in DMQMC, where the distribution along the diagonal approaches zero with β. 55 When H 0 is based on Hartree-Fock, H 0 ii = H ii , the initial condition must also be changed and this is described in detail in the original paper. 55 The grand canonical density matrix corresponding toĤ 0 is used to obtain the desired distribution according to e −βĤ 0 .

C. Symmetric versus asymmetric spawning
For this study in particular, it is important to distinguish between symmetric and asymmetric modes of spawning. In DMQMC as canonically formulated (Eq. (2)) spawning is allowed both on rows and columns because the propagator is symmetric. When we refer to DMQMC in this manuscript, we generally mean this canonical formulation unless otherwise specified. However, it is also possible to have asymmetric DMQMC, with a propagator: where the spawning is restricted on rows (or equivalently on columns). The propagator in IP-DMQMC is canonically asymmetric, with the same spawning restriction as asymmetric DMQMC. While a symmetric propagator exists for the uniform electron gas, 56 it does not for molecular systems and is complicated to develop and test.

D. Full configuration interaction quantum Monte Carlo
Next, we describe the FCIQMC method 29 in brief as it will be used for comparison throughout this work.
FCIQMC begins with the imaginary time Schrödinger equation where |Ψ 0 is the ground state wavefunction,Ĥ is the Hamiltonian operator and τ represents imaginary time.
Here, the wavefunction is represented as a sum over Slater determinants, |D i , where c i is the coefficient on the i th determinant and the Hamiltonian is represented as In the same vein as DMQMC, we can obtain a finite difference equation by substituting the sum over Slater determinants (from Eq. (10)) into Eq. (9), where c m i is the coefficient of the i th determinant at iteration m of the simulation. Note here that the total population of particles, N w , is given by N w = i |c i |. To obtain an estimate of the ground state energy, S is varied to keep the particle population constant, and can be averaged to obtain the estimate.
The rules for evolving particles are those on which the DMQMC algorithm was subsequently based. At each step of the simulation, the particles on each element will undergo spawning, death/cloning and annihilation, as they do in DMQMC. Particles spawn from site i with weight c i to connected sites, j, where i = j, where the probability is uniform in j. In the death/cloning step, particles on site i increase or decrease their population according to |S − H ii |τ . Particles on site i with opposite signs are removed from the simulation.
The particle population is evolved through imaginary time following the rules above, through a system dependent number of iterations. After the wavefunction emerges, the correlation energy is found by averaging over the iterations in the simulation, in a similar fashion to how β loops are averaged in DMQMC.

E. Initiator Adaptation
The sparsity of the thermal density matrix (or the wavefunction coefficient matrix in FCIQMC) can be utilized through use of the initiator approximation variation of both methods, here represented as i-DMQMC 34 and i-FCIQMC 33,57 . The initiator approximation works by setting a threshold "n add " value, where spawning to unoccupied matrix elements is limited to occur only from matrix elements with particle populations larger than n add , called "initiator determinants" (or from coincident spawns of particles of the same sign from two non-initiator sites). This approximation limits the number of density matrix elements (or vector elements in FCIQMC) that need to be sampled over the course of the simulation. Increasing the total number of particles, N w , can reduce the magnitude of the approximation. Both of the original algorithms are obtained as N w → ∞. The initiator adaptation can be used with or without the interaction picture.

F. Kernel Density Estimation
The plateau height in this work is defined as the population that occurs with the highest frequency in the simulation, and we call this population the critical population (N c ). The Scott kernel density estimation (KDE) method 58 is used in this work to assign critical populations through a systematic and reproducible protocol. This is a continuous adaptation from prior work. 59 The KDE method works by calculating the probability that a certain walker population is present in the DMQMC simulation through the use of a KDE kernel, K. If we let f (x) be a continuous function representing the population dynamics in one trajectory, we can use the kernel density estimatorf where h is a smoothing parameter and n is the number of data points to find the KDE kernel. The KDE kernel itself gives a probability distribution of the number of walkers as a function of the number of walkers. The maximum value of the kernel will correspond to the critical walker population. Simulations to measure the plateau height are performed with a single β loop, and the output files are analyzed using the Python scripts provided in the HANDE software package, 51 producing one analysis file per output file. The plateau assignments are performed on the data sets with the total walker population (N w ) on a logarithmic axis. For the plateau heights in DMQMC, there were some cases where the simulations entered variable shift before the annihilation plateau occurred, or the total population collapsed to zero and did not recover. If either of these situations occurred in the simulation, it was not used when measuring plateau heights.
The maximum value of the KDE kernel is assigned as the critical population, and these are collected in a separate file. Graphs of the KDE kernel and the total walker population are produced, and checked visually to ensure the critical population was assigned correctly. Once all plateaus have been validated by visual inspection, the critical populations are averaged and the standard error calculated.
We note here that the FCIQMC critical walker populations used throughout this work are from Ref. 29, and were not recalculated for this work.

III. RESULTS AND DISCUSSION
Calculations were performed on a variety of linear hydrogen chains, and other small atoms and molecules using the HANDE-QMC package, version 1.4 and 1.5. 51 All simulations in this work were performed using a timestep of 0.001 and a shift damping value of 0.30. Integral dump files were generated using MOLPRO, 60 in the form of an FCIDUMP. 61 The single particle eigenvalues for the systems are then calculated using an in-house code from the orbitals in the FCIDUMP according to standard equations. 62 These single particle eigenvalues are then added to the FCIDUMP before the core Hamiltonian energy.
The equilibrium H n chains used in this study had a bond length of 0.945110567 Angstroms, and the stretched H n chains had a bond length of 1.270025398 Angstroms. These correspond to 1.786 a.u. and 2.4 a.u. respectively, which come from a previous study using auxiliary field quantum Monte Carlo. 6 The H 2 O system used a O-H bond length of 0.975512 Angstroms and an H-O-H angle of 110.565 degrees. The CH 4 system used a C-H bond length of 1.087728 Angstroms, and a H-C-H bond angle of 109.47122 degrees. The bond lengths for the diatomic systems are as follows: HF, 0.91622 Angstroms; NaH, 1.885977 Angstroms; C 2 , 1.27273 Angstroms; N 2 , 2.068 a.u. and stretched N 2 , 4.2 a.u. These come from a previous study using FCIQMC. 29 The critical walker populations, or "plateau heights", were measured by the Scott KDE method 58 using NumPy 63 in Python3. The DMQMC calculations to measure critical populations for the H 6 systems were performed with initial populations of 5×10 2 and target populations of 5×10 6 , and for H 8 , the simulations used initial populations of 5 × 10 4 and target populations of 5 × 10 8 . These simulations were propagated to β =25. The IP-DMQMC and FCIQMC simulations for measuring the critical populations were performed with an initial population of 1 and a target population of 5 × 10 8 . All calculations used the integer walker algorithm in all methods to maintain comparability with the plateaus reported in the first FCIQMC paper. 29 . When one walker is used, this means we are sampling exactly one row per β-loop (for asymmetric methods).
The following results are now arranged as follows: in Sec. III A, we begin by confirming the presence of the annihilation plateau and compare the critical walker population in DMQMC to FCIQMC for stretched H 6 , which essentially reproduces known results from Blunt et al. 54 . In Sec. III B, we then explore the connection to the unphysical Hamiltonian related to FCIQMC. 31 Next, we generalize our finding from Sec. III A to a wide range of atomic and molecular systems in Sec. III C, and also explore the interaction picture variant of DMQMC. We then discuss similarities and differences between DMQMC and FCIQMC in Sec. III D, and energy convergence in Sec. III E. Finally, we compare the initiator adaptions to IP-DMQMC and FCIQMC in Sec. III F.
Throughout the manuscript, DMQMC refers to symmetric DMQMC. This is the only type of DMQMC mentioned in Sec. III A, Sec. III B, and Sec. III C. In Sec. III D, asymmetric DMQMC is introduced and discussed and we continue to use it throughout the paper. IP-DMQMC uses asymmetric propagation throughout. In section headings and the captions of figures, information about whether DMQMC is being propagated in a symmetric or an asymmetric fashion is repeated for emphasis and clarity.
A. An example of a symmetric DMQMC annihilation plateau We first begin by describing and then reproducing the original finding of the DMQMC annihilation plateau, where we offer an example of an ab initio system. This section is intended to introduce readers to features of an annihilation plateau. The first paper on DMQMC 54 described the sign problem in this method as similar to that of FCIQMC, due to the close similarities between the population dynamics within the method. Of particular interest is the annihilation step, which is identical between the two methods, and is found to be key in overcoming the sign problem, as described earlier. One difference between the two methods is that the rate of annihilation is likely less frequent in DMQMC than in FCIQMC because there are more density matrix elements than there are terms in the wavefunction coefficient vector. Blunt et al. 54 suggested that because of this slower rate, a higher number of walkers would be needed in DMQMC to overcome the sign problem -approximately the square of the size of the FCIQMC critical walker population -and this observation was based on a Heisenberg model calculation.  The annihilation plateau for stretched H 6 /STO-3G is shown in Fig. 1. This plateau occurs after the first exponential growth, when the population reaches a systemspecific population of walkers, as seen in Fig. 1(a) (for one β loop) between β = 0 and β = 5. When the specific population of particles is reached, the spawning and annihilation rates are approximately equal, resulting in no population growth, i.e., the plateau. After exiting the plateau around β = 15, we observe a second exponential growth phase. The plateau can almost always be visually identified by its distinctive appearance, although in practice, we have also automated this initial measurement (see Sec. II).
When inspecting Fig. 1(b), the instantaneous energy estimate begins the simulation in reasonable agreement with the finite temperature full configuration interaction (ft-FCI) energy, but quickly thereafter the energy fluctuates considerably. After the simulation exits the plateau region, we see a return to agreement between the DMQMC energy and the ft-FCI energy.
We can also compare the plateau heights for stretched H 6 /STO-3G (200 determinants) in DMQMC and FCIQMC. Here, the plateau height is measured at 2.927(5) × 10 4 particles for DMQMC. This system in FCIQMC has a smaller plateau height, at only 2.2(1) × 10 2 particles. This is consistent with the description of Blunt et al., 54 where the authors commented that the DMQMC plateau height is approximately the square of the FCIQMC plateau height. The plateau occurs between β = 5 and β = 15 in this simulation and this temperature range is something that we cannot easily control as an independent variable. Thus, while the critical temperature is something we could measure we generally neglect it for this study.
In summary, the annihilation plateau in DMQMC for ab initio systems follows previous observations based on model Hamiltonians 54 and the FCIQMC annihilation plateau.

B. Connection between the plateau in symmetric DMQMC and the unphysical Hamiltonian and annihilation rate
The sign problem arises in DMQMC because spawning events are affected by the sign of the Hamiltonian, H ik , connecting two density matrix elements. In general, the sign of the matrix element H ik (i = k) can be positive or negative. One way to think about how this arises is that the Slater-Condon rules are applied by bringing the occupied orbitals in determinant i and k into maximum coincidence by permuting the electron indices with each permutation causing a change in sign. It follows that ρ kj can also have any sign. A sufficient number of walkers must be present to allow for the efficient cancellation of signed spawning events arriving at ρ kj to resolve the sign of ρ kj .
Spencer et al. 31 proposed that the sign problem in FCIQMC was: (1) due to an unphysical Hamiltonian (H) whose off-diagonal matrix elements have been wholly negated while leaving the magnitude unchanged, i.e., where δ ik is the Kronecker delta, and (2) tending to be as severe as the energy of the dominant eigenvalue of the unphysical Hamiltonian. The authors found that these can be summarized by the following equation for the critical walker population, N c : Here, κ is the annihilation rate constant and V max is the energy of the highest energy eigenstate of V = −H accounting for the shift correlation energy and the HF energy (i.e. V max = V 0 + S + E HF ). The approximation in the equation refers to this being valid in the limit of a small population and to first order in V max . It will also be helpful to define a variable T max = T 0 + S + E HF . This is the highest energy eigenstate of T = −Ĥ shifted by the same amount as V max . Below, we test the same observations for DMQMC using the stretched H 6 /STO-3G system. No. of Particles (walkers) It is first useful to identify the sign structure of density matrix for both the physical and unphysical Hamiltonians. These are shown in Fig. 2(a) and Fig. 2(b), respectively for β = 3. It can be seen in this figure that these matrices differ in both the signs of their elements as well as the distribution of the occupied elements. In the physical Hamiltonian, there is a mixture of both positively and negatively signed elements, distributed densely across the entire matrix. The combination of the heavily signed and densely packed elements explains why this inverse temperature is difficult to sample. In contrast, we see in the unphysical Hamiltonian matrix that only positively signed elements exist, and is more evenly distributed compared to the physical Hamiltonian matrix. Now, in the DMQMC simulations of both Hamiltonians, different dynamics are seen. For the physical Hamiltonian, DMQMC exhibits a characteristic plateau shape as the total population growth rises exponentially, pauses, and then resumes (Fig. 2(c)). Only when the population growth resumes does the growth of walkers on the diagonal of the density matrix start in earnest. In the dynamics of the simulation, we see that the walkers on the diagonal tend to spawn and then die, depleting the diagonal population. It is only when enough of a population exists on the off-diagonal part of the density matrix and the sign structure has been established that the diagonal population can be sustained. By contrast, for the unphysical Hamiltonian, DMQMC exhibits largely uninterrupted growth in both the total population and the population of walkers on the diagonal. In this case the sign problem has occurred becauseH =Ĥ; this condition is necessary but not sufficient. However, there are examples where it is sufficient to have a similarity transformation which maps the two matrices onto each other. Such is the case in the bipartite Heisenberg model, for example. 31 To show that the sign problem is also related to the dominant eigenvector of the unphysical Hamiltonian, ft-FCI results are shown in Fig. 2(d). We can see here that in general, the energies obtained from the two Hamiltonians are different, where the unphysical Hamiltonian energy is lower than that of the physical Hamiltonian energy. The one exception we see is at β = 0, when the two solutions are degenerate owing to the trace being the same between the physical and unphysical Hamiltonian. In the low β regime is exactly where the dynamics appear to be the most similar in terms of population growth (Fig. 2(c)) which is consistent with the idea that the dominant eigenvalue of the unphysical Hamiltonian causes a change in the population dynamics.
To analyze this further, we can compare the population growth rates when usingĤ andH. Assuming a growth rate of N w ∼ e kβ , we can find the instantaneous rate constant for growth from d dβ ln(N w ). This is shown in Fig. 3. The growth rate for theH propagator oscillates around e Vmax β for the whole of the simulation. By contrast, the growth rate forĤ in the pre-plateau region starts at e Vmax β, while post-plateau the growth rate tends towards e Tmax β at large β. This lends further evidence to the relationship between the pre-plateau dynamics and H.
To provide further data to make the point that the dominant eigenvalue causes a change in dynamics, we scaled the off-diagonal matrix elements linearly by a factor of C, starting from the true Hamiltonian, resulting inH where H ik and δ ik follow previous definitions, with an additional factor of a positive constant C. The plateau height for low C follows a linear trend (Fig. 4) which fits the form of Eq. (14) as V max is linear in C for this system (assuming a constant κ). This observation is also consistent with that of Spencer et al. 31 showing that the plateau height varies linearly with U/t in the Hubbard model where U is the on-site interaction strength and  t is the hopping integral. Thus the analog to U in our re-scaled molecular Hamiltonian is C.
For completeness, the last component of the plateau expression given in Eq. (14) we want to test for DMQMC is the dependence on the shift parameter, S. We collected data for the equilibrium H 8 system shown in Fig. 5. It can be seen from these data that at low S, the plateau height is linear in S which is consistent with the form of Eq. (14).
In this section, we found that the sign problem and population dynamics in DMQMC can be related to similar observations made of FCIQMC. 31 In the next section we explore the relationship between the plateau heights of the two methods along with IP-DMQMC.

C. How the symmetric DMQMC and IP-DMQMC plateau heights scale in relation to the FCIQMC plateau height
In Sec. III A, we observed that the plateau height in DMQMC was approximately the square of the plateau height in FCIQMC for the stretched H 6 system. In this section, we attempt to generalize this observation to a wide range of atomic and molecular systems for both DMQMC and IP-DMQMC (which was outlined in Sec. II B). To achieve this, we study the range of closed-shell systems 67 that were previously considered by Booth et al. 29 (various atoms and molecules comprised of first-row atoms) and supplement these with 1D hydrogen chains. The latter set are of interest because they are approximate analogs of the Hubbard models, which are also used for plateau studies 31,59 . Thus, the total test set is comprised of: Ne (aug-cc-pVDZ), H 2 O (cc-pCVDZ), HF (cc-pCVDZ), NaH (cc-pCVDZ), C 2 (cc-pVDZ), CH 4 (cc-pVDZ), N 2 (cc-pVDZ), stretched N 2 (cc-pVDZ), as well as stretched and equilibrium H n (STO-3G) for even n between 4 and 16 inclusive. The Be atom is excluded from the test set as it has no measurable annihilation plateau in IP-DMQMC. This test set represents a variety of chemical systems including hetero-and homonuclear diatomics with single and multiple bonds. Our preliminary observation was that the DMQMC and IP-DMQMC plateau heights were a system-dependent fraction of the size of the space similar to FCIQMC. This made it difficult to establish a specific trend with system size.
We anticipate that each system will have a plateau height which is a system-dependent fraction of the size of the space, similar to FCIQMC. We therefore plot the DMQMC plateau height against the FCIQMC plateau  29 (circle), equilibrium Hn chains (even n between 6 and 16 inclusive, square symbol) and stretched Hn chains (even n between 6 and 16 inclusive, × symbol). IP-DMQMC simulations used a target β = 25. Straight lines are plotted for both y = x (solid) and y = x 2 (dashed) to help guide the eye. The plateau heights were measured using the KDE method, and were averaged over 25 simulations. The FCIQMC critical walker populations are from published data. 29 Error bars are shown, and, in some cases, are smaller than the size of the marker. In this figure, DMQMC is symmetrically propagated and IP-DMQMC is asymmetrically propagated. height for the same system ( Figure 6). These values were available for equilibrium and stretched H 6 , and for equilibrium and stretched H 8 . All of the other systems in our test set had critical populations in DMQMC that were > 5 × 10 8 particles (our choice of the cutoff in population in our experimental design). What we see in this data is that for these four systems, the DMQMC plateau height is approximately the square of the FCIQMC plateau height.
We now turn our attention to the interaction picture variant of DMQMC (IP-DMQMC). While this was introduced in Sec. II B, it is instructive to provide a number of details at this point. IP-DMQMC targets a specific β value (here β = 25 to consistently allow the plateau to be found) and initializes on an exactly known auxilliary matrix, (f (τ ) = e −(β−τ )Ĥ 0 e −τĤ ), with the weights of the auxilliary matrix replacing the random sampling of the diagonal identity matrix in DMQMC. IP-DMQMC also modifies the propagator such thatf (τ = β) = ρ(β) and that the propagation is asymmetric (i.e. only happens down the rows or columns of the density matrix). Figure 6 also shows plateaus heights from IP-DMQMC. No. of States N c (IP-DMQMC) N dets Figure 7. The critical population (Nc, blue) in IP-DMQMC for the atoms and molecules in our test set, compared to the estimated size of space (N dets , red). These Nc were collected in the same way as Fig. 6. In this figure, IP-DMQMC is asymmetrically propagated.
Our calculations in Fig. 6 show that the IP-DMQMC plateau height is approximately equal to the plateau height in FCIQMC for the systems studied here. For example, for the stretched H 6 system, the IP-DMQMC plateau height is 2.2(1) × 10 2 particles, and the FCIQMC plateau height is also 2.2(1) × 10 2 particles. This finding is remarkable, as it shows that the critical walker population in IP-DMQMC is directly related to the same in FCIQMC.
To further emphasize this, Fig. 7 shows the critical walker population in IP-DMQMC plotted next to the size of the Slater determinant space in FCIQMC. It can be seen that almost all of these systems have plateaus heights lower than the number of determinants and, therefore, lower than the square root of the number of elements in the density matrix.

D. How IP-DMQMC has the same plateau height as FCIQMC
In order to examine what differences in the DMQMC and IP-DMQMC methods give rise to different critical populations, we begin by analyzing Eq. (14). If we assume that V max is the same or approximately the same then κ can be calculated for each method. For the sim- ulations of stretched H 6 , we can find V max = 1.677 Ha by diagonalization. Using the plateau height, we can then find that the κ values for FCIQMC, DMQMC, and IP-DMQMC are 7.3×10 −3 , 5.7×10 −5 , and 7.3×10 −3 respectively. Here, we can see that the IP-DMQMC rate of annihilation is the same as FCIQMC and approximately the square root of that in DMQMC i.e. IP-DMQMC requires a similar rate of annihilation as in FCIQMC to resolve the sign problem. These can be corroborated by measuring the large-β limit growth rate of the population in Fig. 2(c) and through measuring the annihilation rate directly from the number of walkers removed in the simulation. The walkers removed by annihilation are shown in Fig. 8. The graph shows agreement with the observation above, that the annihilation rate agrees between FCIQMC and IP-DMQMC, and both are much lower than the rate in DMQMC. Going a step further, we can show that IP-DMQMC and FCIQMC have more similarities. Most notably, when IP-DMQMC is started from one walker, the propagator reduces to that of FCIQMC, exactly. To demonstrate this, we start with the IP-DMQMC propagator from Sec. II, recalling thatĤ 0 is diagonal. If we assume that our one walker lands on the zeroth row, thenfĤ 0 = H 00 f 00 and will only affect the diagonal. Then, the contribution to ∆f which is equal toĤf leads tof = (1 +Ĥ)f , which is the FCIQMC propagator. The element H 00 refers to D 0 |Ĥ|D 0 = E HF . In IP-DMQMC the term H 00 f 00 modifies the Hamiltonian, subtracting the Hartree-Fock energy from the propagator, as in FCIQMC. This particular similarity between IP-DMQMC and FCIQMC is what guarantees the The critical population of different rows in the stretched H6 density matrix as a function of their Vmax value. Each row has its own critical population, and was measured from one β loop. The zeroth row is marked with a red × symbol. Simulations for measuring the plateau height were started with one walker and had a target population of 5×10 8 . Vmax was calculated using an in-house analytical IP-DMQMC code, by propagatingH to β = 25 separately for each row. The energy was calculated at the beginning and end of the simulation, and Vmax is equal to the difference between the final and initial energies. The linear fit (green dashed line) is given by Nw(β) = 8.2(4) × 10 2 (Vmax) − 1.4(1) × 10 3 . In this figure, IP-DMQMC is asymmetrically propagated.
equivalence of the critical populations in Fig. 6, provided that the zeroth row off (in IP-DMQMC) is only chosen during initialization. It is reasonable to assume that when a high target β value is used (as in our simulations shown in Fig. 6), the zeroth row will indeed be chosen. Thus, Fig. 6 only represents a minimal plateau in IP-DMQMC when the ground-state outer product is being simulated. Unless the simulation is run at very high β, we can expect that other rows will need to be simulated. Other rows are not encountered during an IP-DMQMC simulation started from one walker because the propagator prevents other rows from being accessed during the simulation. This means that we can also measure N c on a per row basis. When IP-DMQMC is deliberately initialized on different rows, we find that there are slight changes in the plateau as we move away from the zeroth row. In order to understand these changes, we note that N c ∝ V max but that the effective V max for a given row requires the −Ĥ (0)f term in the propagator is taken into account. In practice, this means that the effective V max is raised by |H ii | for row i. The critical populations from different rows in IP-DMQMC are shown in Fig. 9 for stretched H 6 , showing the linear relationship predicted by N c ∝ V max . The average critical population (taken as an average over rows) is N c = 1.22(3) × 10 3 , which is slightly higher than reported in the previous section (N c = 2.2(1) × 10 2 ). The larger plateau height is due to the influence of the −Ĥ (0)f term in the IP-DMQMC propagator raising the plateau relative toĤρ the unmodified propagator (used in DMQMC).
One last question we consider in this section is the The total walker populations (Nw(β)) for IP-DMQMC (blue) and asymmetric DMQMC (green) as a function of inverse temperature (β) for a random row that is not the zeroth row from Fig. 9. In the asymmetric DMQMC simulation, the shift was set to Hii to match the IP-DMQMC methodology. The critical populations for these IP-DMQMC and asymmetric DMQMC simulations are shown as markers for the two methods (black × symbol and black triangle, respectively). Simulations were started with one walker. In this figure, IP-DMQMC and DMQMC are both asymmetrically propagated.
following: Does the interaction picture or the asymmetric propagation cause IP-DMQMC to have FCIQMClike plateau heights? Until now, we have been running DMQMC in symmetric mode, where rows do interact because spawning occurs along rows and columns. Recall that earlier in this section, we showed that IP-DMQMC has a lower annihilation requirement than DMQMC. In practice, we find when DMQMC is run in asymmetric mode, it mirrors IP-DMQMC in having a lower plateau height than symmetric DMQMC. We actually find that the two are identical, if the diagonal shift in asymmetric DMQMC matches what is subtracted off byfĤ 0 in IP-DMQMC. This is shown by visual inspection in Fig. 10. We tested the difference by running 25 β loops. We found the average plateau heights to be N c = 3.7(2)×10 2 and N c = 3.44(9) × 10 2 for IP-DMQMC and asymmetric DMQMC respectively and the difference is not statistically significant. We believe, therefore, that the choice of spawning mode (symmetric versus asymmetric) gives rise to the differences we see between different methods in Fig. 6.
In this section, we have determined the differences between DMQMC and IP-DMQMC that gives rise to different critical populations. A key finding was that, because the critical population measurement is started from one walker, IP-DMQMC and asymmetric DMQMC only has one row occupied throughout the whole simulation (due to the structure of the propagator). For large target β this is likely to be the ground-state-like row mak-ing one-row IP-DMQMC equivalent to FCIQMC. For the non-ground-state rows, each can have a shift applied to make the plateau equivalent but, without modification, the plateau grows slightly. A one row asymmetric simulation on its own does not allow for a reliable thermal energy to be obtained. To obtain an accurate thermal energy, an average over row simulations must be found e.g. by using β loops. The number of β loops required to converge the energy thus plays a role in the scaling of IP-DMQMC beyond the sign problem. Overall, then, this has the effect of allowing the distribution of memory costs across different β loops, perhaps allowing for the convergence of systems that are too large for symmetric propagation. However, the question remains as to how efficient this averaging is and whether there is a gain in cost relative to the DMQMC plateau. The subject of the stochastic error encountered when sampling the different rows of the density matrix is the subject of the next section.

E. Energy convergence with respect to number of rows and β loops in IP-DMQMC and asymmetric DMQMC
We now attempt to work out how many rows are required to converge an asymmetric DMQMC or IP-DMQMC calculation. It is, in principle, possible to converge a calculation either using row sampling (from the starting point of the simulation) or more beta loops. We use an analytical implementation of IP-DMQMC and asymmetric DMQMC to measure the energy convergence of stretched H 6 with respect to N rows × N β by carefully controlling the type of sampling, the number of rows and the number of β loops. In the analytical code, walkers are initialized by randomly placing walkers (one at a time) on diagonal elements of the density matrix. The random distribution is uniform for asymmetric DMQMC and normalized thermal weights are used for IP-DMQMC. The propagation steps are handled deterministically, which removes the sign problem and allows us to isolate how many rows need to be sampled. The error is calculated in the normal way, using analysis tools provided in the HANDE-QMC package.
In general, the energy was well converged within error bars across the whole of the data set. This is in part because the H 6 system contains only 200 rows (or that there are 200 FCI determinants) which means we were oversampling in general. However, even for N rows ×N β < 200, we see that the energy is converged within error (< 2σ) for the majority of the data set 68 . One example of this is when a minimal number of rows is sampled, which shows that the rows can be sampled independently in IP-DMQMC and asymmetric DMQMC (Fig. 11(a)). This is important because it at least means the memory requirement of the plateau storage (lowered due to asymmetric propagation) can be distributed across different instances of IP-DMQMC or asymmetric DMQMC as suggested in the previous section.
In general, we found that there was an trade-off between N rows and N β when it came to reducing the stochastic/sampling error. This can be seen in graphs of the stochastic error plotted against N rows × N β , where all of the data sets are (by visual inspection) part of the same distribution. This distribution generally decays according to a power-law fit of (N rows × N β ) in the large N rows or N β limit. Examples of this are shown for two representative β values in Fig. 11(b) and Fig. 11(c). On these graphs, the error has been multiplied by (N rows ) to make fair comparison between IP-DMQMC and asymmetric DMQMC. In the graphs shown we also see that β = 1.0 generally has a higher error than β = 7.0. For β = 7.0, IP-DMQMC has lower stochastic error than asymmetric DMQMC by an order of magnitude while at β = 1.0 their error is more comparable. This advantage appears to be reduced at low N rows × N β , which is the limit we want to be able to run our calculations in. It is still possible to see (in Fig. 11(a)) that IP-DMQMC has a lower systematic error, indicating that the asymmetric DMQMC error may have an under-sampling error.
Overall, we find that it is possible to take maximal advantage of the reduced plateau height of IP-DMQMC by running simulations on individual rows of the density matrix and averaging over β loops. For higher temperatures (or for asymmetric DMQMC), the whole diagonal of the density matrix must be sampled which is likely to be costly. However, for lower temperatures (higher β) IP-DMQMC generally converges with a smaller systematic difference to the exact result, and a smaller stochastic error. IP-DMQMC will exhibit reduced computational cost compared to asymmetric DMQMC due to a need to sample fewer rows (whether through walkers or β loops).

F. The initiator approach applied to IP-DMQMC
The initiator adaptation (Sec. II E) was developed to maintain population on the diagonal which is popular in FCIQMC because it removes the requirement that the simulation has to have a total walker number greater than the critical walker population (i.e. the plateau is removed), introducing only a modest error. Unfortunately, the removal of the plateau means we cannot compare how i-FCIQMC and i-IP-DMQMC scale using this measure alone. We can, therefore, use a previous study of i-FCIQMC 57 where a walker population threshold measure (N thresh ) of 50,000 walkers on the Hartree-Fock determinant was used as the population requirement for  a converged simulation. This threshold is analogous to measuring the plateau height because it was shown for a variety of atoms that the energy did not vary after this threshold was reached and the simulation was converged with respect to stochastic sampling, 57 which is the same idea as the canonical method needing to reach a critical walker population to obtain a converged energy. We considered but did not attempt a "growth witness" measure (G) introduced by other authors. 64 To adapt this measure for IP-DMQMC, we consider a threshold of 50,000 walkers on the trace of the density matrix which controls for systematic and stochastic errors simultaneously. We found that this threshold provides simulations with a mostly consistent stochastic error across system sizes ( Fig. 12(a)). In this section, we compare i-IP-DMQMC with i-FCIQMC.
In Fig. 12, we show the total walker population at the simulation iteration at which the population threshold was met on the diagonal of the matrix (or HF determinant for i-FCIQMC), and we will refer to this walker value as N thresh throughout this section. Here, we find that the i-IP-DMQMC cost in terms of walker number is the same as i-FCIQMC for low temperatures (i.e., β=10), which makes sense because we are simulating the ground state when our choice of target β is large. In the intermediate temperature range (β=2 and β=5), we find that for the smaller hydrogen chains, the cost is slightly higher in i-IP-DMQMC, but as the length of the chain is increased, i-IP-DMQMC returns to being approximately the same cost as i-FCIQMC. At the two higher temperatures (β=1 and β=0.1), we find the cost is lower than that of i-FCIQMC. We attribute this to the initiator adaption itself. This variation tends to keep walkers on the diagonal of the density matrix for IP-DMQMC and at high temperatures more particles on the diagonal is closer to the physical solution, making it easier for IP-DMQMC to simulate. In general, these data show that the initiator approximation in IP-DMQMC controls the walker population in a similar manner to the initiator approximation in FCIQMC. This gives us confidence that the initiator approximation can be used in future applications.
We note in passing that the importance sampling of DMQMC was also designed to maintain particles on the diagonal of the density matrix, 54 and is explored in Appendix B. We also note that we have only looked at stochastic error here and a study of systematic initiator error is extremely important going forward. Due to its complexities in the DMQMC method, 56 this is left for a future study.

IV. CONCLUSIONS
DMQMC has been shown to be a promising method for finite temperature applications, and in this work we have confirmed that DMQMC (especially in its interaction picture variant) shows the potential to be as effective  Figure 12.
(a) The stochastic error in i-FCIQMC and i-IP-DMQMC for linear Hn chains in the STO-3G basis set, for n = 4, 6, 8, 10, 12. Simulations of both methods were performed at increasing target populations starting at 100 walkers and increasing to 5 × 10 6 walkers. For each method, a single simulation was used to determine the smallest target population (N thresh ) required to reach 50,000 walkers on HF. For i-IP-DMQMC, the stochastic error shown was obtained by averaging over 5 β loops that reached N thresh . For i-FCIQMC, the stochastic error was found for 5 different simulations (50,000 report cycles, a timestep of 0.001 and an initial population of 10 walkers), and then averaged. Error bars show one standard error. No shift damping was used in these simulations. (b) The total walker population at the simulation iteration at which the population threshold was met on the diagonal of the matrix (HF determinant for i-FCIQMC) in the i-IP-DMQMC simulation (N thresh ) as a function of the same in the i-FCIQMC simulation (N thresh ). with initial populations of 10 particles, and with target β values of 0.1 (magenta × symbols), 1 (blue circles), 2 (green diamonds), 5 (gold squares) and 10 (purple + symbols). The y = 10 b × x m fits are shown in the legend, on the same line as the marker corresponding to the β value. In this figure, IP-DMQMC is asymmetrically propagated.
for finite temperature work as FCIQMC is for ground state simulations. We confirmed that the critical walker population in symmetric DMQMC scales as the square of that in FCIQMC. Additionally, we found that the critical walker population in IP-DMQMC is the same as that of FCIQMC across all β values due to the asymmetric sampling in IP-DMQMC. We also determined that the trade-off between sampling a small amount of rows many times versus sampling all rows fewer times is approximately equal, opening an additional avenue of development for the method. The latter is a very exciting result, as it shows that we can obtain a temperature-dependent energy at roughly the same memory and walker cost as FCIQMC, allowing us to treat systems with IP-DMQMC that cannot be treated by DMQMC. With respect to the critical walker population, this implies that IP-DMQMC has more utility compared to DMQMC, because a smaller population of particles is required to obtain the density matrix associated with the physical Hamiltonian. Finally, we showed that the initiator adaption with IP-DMQMC performs in a similar way than i-FCIQMC, again allowing us to expand upon the systems we can treat with this method.
As such, we now know that IP-DMQMC will be more useful than DMQMC for systems with a severe sign problem. One disadvantage of using IP-DMQMC, which we did not explore here, is that IP-DMQMC requires separate simulations to obtain energies for different inverse temperature values. Whether the computational overhead is then more expensive to obtain a full β spectrum in IP-DMQMC compared to DMQMC is still an open question. We note in passing that we did not explore the connection between this observation and the Krylov projected FCIQMC, 65 as we felt that this was beyond the scope of this work.
Overall, this strongly suggests that the IP-DMQMC algorithm has the same potential as FCIQMC, and gives a focus for future development. A natural place for future work to begin is to explore the uses (and limitations) of the initiator approach in a systematic way as well as examining ways to modify and lower the IP-DMQMC plateau height. For example, as it is known that basis function rotations do affect the plateau height, we are inclined to explore basis functions that are optimized for a given temperature. 7 IP-DMQMC differs from DMQMC in that a specific (target) β must be specified and then each β has a unique simulation. This means that there is the potential for dependence of the IP-DMQMC plateau height on the specified β. In this section, we explore whether the IP-DMQMC plateau heights across intermediate β values are the same as the FCIQMC plateau height.

V. ACKNOWLEDGEMENTS
We test this dependence using the H 2 O system, and we present critical populations as a function of target β in Fig. 13. As each target β progresses through the plateau to a varying extent by the time the target is reached, to ensure a fair test we used a wall time limit of 4 hours instead. We found that for all β values simulated the plateau heights are between 4 × 10 7 and 5 × 10 7 , showing evidence that the IP-DMQMC critical population is not strongly β dependent. The FCIQMC critical population for this system was found to be 4.74 × 10 7 , and so these results confirm that the IP-DMQMC critical population is approximately the same as FCIQMC across all β values.
When moving to a different system, namely stretched H 6 , we found that it was more challenging to measure a plateau height at intermediate β values directly, as changes in input parameters would be required to make sure that a plateau even exists by the time the target β is reached. When comparing the two systems, the walker growth is slower in the simulations of the stretched H 6 system compared to the growth in the simulations of H 2 O, and we expect this may be why the plateau heights in H 2 O can measured directly, whereas they cannot be measured directly in stretched H 6 . By means of an alternative, we instead study how the (symmetric) DMQMC energy converges to the exact result with walker number and how this convergence changes a function of the target β. We do so as an alternative to changing input parameters in the IP-DMQMC simulations.
To test these values, simulations were performed at varying populations between 10 2 and 10 6 walkers with variable shift used throughout the simulation. We note that the random initialization algorithm of IP-DMQMC means the population can vary slightly from the starting population. For this test, the number of β loops was reduced as the walker number was increased so the error stated is then the stochastic error for a given computational cost. By visual inspection of the energy differences to ft-FCI (Fig. 14), the energy differences rapidly fall to zero above N w = 10 5 walkers for DMQMC and N w = 10 3 walkers for IP-DMQMC. After this, energies are well converged with relatively small error bars. This is consistent with the plateaus estimated at large β referenced above. In the case of DMQMC below N w = 10 5 walkers, the noise grows with rising β, which is consistent with an exponentially falling signal-to-noise ratio which characterizes the sign problem. In particular, for DMQMC simulations below the plateau, the trace becomes very small and the energies become difficult to converge.  Figure 14. For the stretched H6 system, the absolute energy difference from exact diagonalization (ft-FCI) is shown for (a) DMQMC and (b) IP-DMQMC. The closer to zero the energy difference is, the more converged we consider the energy. Error bars are shown but may be smaller than the size of the marker in some cases. In both figures, five populations are shown: 10 2 (orange circle), 10 3 (blue star), 10 4 (teal diamond), 10 5 (green × symbol), and 10 6 (purple cross). For DMQMC at Nw = 10 2 , the energy differences at β = 7 and β = 8, are outside the range of the plot, and are −1.814 Ha and −9.808 Ha, respectively. In this figure, DMQMC is symmetrically propagated and IP-DMQMC is asymmetrically propagated.
proach does keep walkers on the diagonal of the density matrix, it also generally slightly raises the height of the annihilation plateau (Fig. 15). We note that this agrees with another study on FCIQMC in the literature. 66 So, while it has promise in terms of converging the energy with reduced noise (due to having an increased trace pop- ulation) we do not investigate it any further here.