The Grand Canonical Ensemble

Our general approach to the quantum many-body problem is to discretize an approximation of the partition function \(Z = \text{tr}[ \exp(-\beta(H - \mu N - h \cdot S))]\), a grand canonical ensemble, via a Suzuki-Trotter decomposition or simply a Trotterization.

For a discretization of space into \(N_x\) sites on a side we regroup all the dimensionful parameters into

\[\begin{split}\begin{align} \tilde{H} &= HML^2 & \tilde{\beta} &= \beta/ML^2 \\ \tilde{\mu} &= \mu ML^2 & \tilde{\vec{h}} &= \vec{h} ML^2 \end{align}\end{split}\]

With an on-site interaction, holding \(\tilde{a}=2\pi a/L\) fixed yields a dimensionless Hamiltonian that depends only on \(N_x\) and the scattering through the AnalyticTuning. This convention makes it simple to understand the Hamiltonian limit (\(N_t\rightarrow\infty\) holding \(\tilde{\beta}\) fixed) in which the Trotterization disappears, and the spatial continuum limit (\(N_x\rightarrow\infty\) holding the physical parameters, and therefore, \(\tilde{a}\), \(\tilde{\beta}\), \(\tilde{\mu}\), \(\tilde{\vec{h}}\) fixed).

The Fermionic Discretization

Trotterizing the thermal circle yields a euclidean time, which can be combined with the Lattice into a SpaceTime.

tdg.Spacetime

class tdg.spacetime.Spacetime(nt, lattice)[source]

Bases: H5able

nt: The number of timeslices.

Lattice: The spatial lattice.

t: The coordinates in the time direction.

sites: The total spacetime (integer) volume.

dims: The spacetime dimensions.

coordinates: A tensor of size [sites, len(dims)]. Each row contains a set of coordinates. Time is slowest and each timeslice matches Lattice.coordinates in order.

vector(*dims)[source]

Parameters: dims (tuple) – Specifies how how many spacetime vectors to produce.
Returns: A dims-dimensional stack of zero vectors. Each vector is of .shape == [nt, Lattice.vector().shape].
Return type: torch.tensor

tdg.FermionMatrix

Fermions are tricky to handle numerically, as they are anticommuting. Instead we introduce an auxiliary field to linearize the action in fermion operators, and then integrate the fermions out, yielding a FermionMatrix. Including the determinant of the FermionMatrix in the sampling includes the fermions, even though the auxiliary fields are bosonic.

Note

The arguments are dimensionless, as above! Henceforth with only write the tilde if clarity requires.

class tdg.FermionMatrix(spacetime, beta, mu=tensor(0.), h=tensor([0., 0., 0.]))[source]

Bases: H5able

The fermion matrix \(\mathbb{d}\) which corresponds to

\[\begin{split}\begin{align} \mathbb{d} &= \begin{pmatrix} -\mathbb{1} & 0 & 0 & 0 & \cdots & 0 & -\mathbb{B}^{-1}\mathbb{F}_{N_t}( h ) \\ +\mathbb{B}^{-1}\mathbb{F}_{1}( h ) & -\mathbb{1} & 0 & 0 & \cdots & 0 & 0 \\ 0 & +\mathbb{B}^{-1}\mathbb{F}_{2}( h ) & -\mathbb{1} & 0 & \cdots i & 0 & 0 \\ 0 & 0 & +\mathbb{B}^{-1}\mathbb{F}_{3}( h ) & -\mathbb{1} & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & -\mathbb{1} & 0 \\ 0 & 0 & 0 & 0 & \cdots & +\mathbb{B}^{-1}\mathbb{F}_{N_t-1}( h ) & -\mathbb{1} \end{pmatrix} \end{align}\end{split}\]

where the double-struck quantities

\[\begin{align} \mathbb{F}_t &= \left(e^{A_{t}+\Delta t\mu + \frac{1}{2}\Delta t\vec{h}\cdot\vec{\sigma}}\right) & \mathbb{B} &= (B = \exp \Delta t \kappa) \otimes \mathbb{1} \end{align}\]

are space ⊗ spin in size.

In the language of Ref. [8] this is called the \(\alpha=0\) exponential discretization because the hopping matrix \(\kappa\) which encodes the kinetic piece of the Hamiltonian appears inside an exponential, and we’re in the attractive channel, so the auxiliary field appears without an explicit \(i\) in \(\mathbb{F}\).

Parameters

spacetime (Spacetime) –
beta (torch.tensor scalar) – The inverse temperature
mu (torch.tensor scalar) – The chemical potential
h (torch.tensor triplet) – The external field

Todo

We should add a method to do fermion matrix-vector multiply, perhaps as __call__. The programming puzzle is that when \(\vec{h}=0\) or \(\vec{h}\parallel \hat{z}\) we can split \(\mathbb{d}\) into two independent pieces that each operate on spacetime vectors, while for generic \(\vec{h}\) we need one matrix that operates on a spacetime ⊗ spin vector.

U(A)[source]

\[U(A) = B^{-1} F_{N_t} B^{-1} F_{N_t-1} \cdots B^{-1} F_{2} B^{-1} F_{1}\]

where \(F_t = \exp A_t\) is a digonal matrix and does not include the chemical potential \(\mu\).

\(B\), \(F\) and \(U\) are space × space, with rather than (space ⊗ spin) × (space ⊗ spin).

Parameters: A (torch.tensor) – An auxiliary-field configuration
Returns: \(U(A)\)
Return type: torch.tensor

UU_tensor(A)[source]

\[\mathbb{U} = \mathbb{B}^{-1} \mathbb{F}_{N_t} \cdots \mathbb{B}^{-1} \mathbb{F}_{2} \mathbb{B}^{-1} \mathbb{F}_{1}\]

where \(\mathbb{F}\) contains the chemical potential and external field terms. However, since those terms are space-independent we can commute them with \(\mathbb{B}\) so that

\[\begin{split}\mathbb{U} = \exp\left(\frac{1}{2}\beta h\cdot\sigma\right) \begin{pmatrix} zU & 0 \\ 0 & zU\end{pmatrix}\end{split}\]

where \(U\) is given by U() and is (space × space).

Parameters: A (torch.tensor) – An auxiliary-field configuration.
Returns: \(\mathbb{U}\) not as a matrix but with shape [space, space, spin, spin]. This makes the spin indices broadcastable, which is often useful for future manipulation.
Return type: torch.tensor

UU(A)[source]

Parameters: A (torch.tensor) – An auxiliary-field configuration
Returns: The same as UU_tensor() but as a true matrix of shape (space ⊗ spin) × (space ⊗ spin), with spin fastest.
Return type: torch.tensor

Todo

Leveraging this makes some code much slower. For example,

>>> UU = torch.stack(
...     tuple(Action.FermionMatrix.UU_tensor(cfg) for cfg in configurations)
...     ).permute(0,1,3,2,4).reshape(-1,
...                                  2*Action.FermionMatrix.Spacetime.Lattice.sites,
...                                  2*Action.FermionMatrix.Spacetime.Lattice.sites)

is much faster to subsequently manipulate than

>>> UU = torch.stack(
...     tuple(Action.FermionMatrix.UU(cfg) for cfg in configuration))

which is something that ought to be understood / addressed. Perhaps it is autograd- or computational-graph-related?

Is it possible to construct this code to distribute over a whole ensemble of auxiliary fields? This perhaps should inform the way we code other observables.

logdet(A)[source]

The log determinant of the fermion matrix appears in the action.

Parameters: A (torch.tensor) – An auxiliary-field configuration.
Returns: \(\log \det \mathbb{d}(A)\)
Return type: torch.tensor

Note

On implementation:

One may show that \(\det \mathbb{d} = \det \mathbb{1} + \mathbb{U}\) where all of the space-independent pieces of \(\mathbb{F}\) can be pulled through so that

\[\begin{split}\mathbb{U} = \exp\left( \frac{1}{2} \beta \vec{h}\cdot\vec{\sigma} \right) \begin{pmatrix} zU & 0 \\ 0 & zU \end{pmatrix}\end{split}\]

Furthermore, an orthogonal transformation can diagonalize the spin-dependent factor and leave the spin-independent factor alone. The determinant of \(\mathbb{d}\) can then be expressed as

\[\det \mathbb{d} = \det\left( 1 + e^{+\frac{1}{2}\beta \sqrt{h^2}} z U\right) \det\left( 1 + e^{-\frac{1}{2}\beta \sqrt{h^2}} z U\right)\]

where the fugacity \(z = \exp(\beta \mu)\) and the \(\sqrt{h^2}\) helps correctly handle complex \(h\).

tdg.Action

The Action is combination of the fermion determinant and the gaussian piece of the auxiliary field.

class tdg.action.Action(spacetime, potential, beta, mu=tensor(0.), h=tensor([0., 0., 0.]), fermion=<class 'tdg.fermionMatrix.FermionMatrix'>)[source]

Bases: H5able

Parameters

spacetime (tdg.Spacetime) – the euclidean spacetime
potential (tdg.Potential) – the interaction with which the fermions interact; must be negative-definite for the Hubbard-Stratanovich transformation used to make sense.
beta (torch.tensor scalar) – the inverse temperature
mu (torch.tensor scalar) – the chemical potential
h (torch.tensor) – the spin chemical potential; a triplet
fermion (tdg.FermionMatrix) – the fermion matrix corresponding to the desired discretization

An auxiliary-field spacetime action for fermions interacting via a potential with inverse temperature \(\beta\), chemical potential \(\mu\), spin chemical potential \(\vec{h}\).

The action is used for importance sampling, since the partition function Z is

\[\begin{align} Z &= \int DA\; e^{-S} & S &= \frac{1}{2} \sum_t A_t (-\Delta t V)^{-1} A - \log \det \mathbb{d} + \frac{N_t}{2} \text{tr} \log \left(-2\pi \Delta t V\right) \end{align}\]

where \(\mathbb{d}\) is the fermion matrix, the time step \(dt = \beta/N_t\), and everything should be understood to be dimensionless.

The last term is a constant normalization to make \(Z\) truly equal to the Trotterization of \(\text{tr} e^{-\beta H}\).

Spacetime: The Spacetime on which the action is formulated.

Potential: The Potential \(V\) with which the fermions interact.

beta: The dimensionless inverse temperature \(\tilde{\beta} = \beta/ML^2\).

dt: The temporal discretization \(dt = \texttt{beta} / N_t\).

mu: The chemical potential \(\tilde{\mu} = \mu ML^2\).

h: The spin chemical potential \(\tilde{\vec{h}} = \vec{h} ML^2\).

V: The spatial representation of Potential on the Spacetime.Lattice

Vinverse: The inverse of Potential on the Spacetime.Lattice

FermionMatrix: The fermion matrix that gives the discretization.

Note

On the chemical potential:

Recall that requiring the contact interaction be written as the quadratic \(nVn\) induces a term in the Hamiltonian proportional to \(n\) itself, which looks just like a chemical potential term. This term comes with a coefficient equal to \(-N_x^2 C_0/2\).

Our sign convention is that the internal energy is \(H-\mu N - h\cdot S\) so the the signs conspire. The fermion matrix is constructed with the ‘offset’ chemical potential \(\mu + N_x^2 C_0/2\).

normalizing_offset: The A-independent contribution to the action needed to match Z to the Trotterized operator definition.

\[\frac{N_t}{2} \text{tr} \log\left(-2\pi \Delta t V\right)\]

gauge(A)[source]

Parameters: A (torch.tensor) – an auxiliary field configuration
Returns: \(\frac{1}{2} \sum_t A_t (-\Delta t V)^{-1} A_t\)
Return type: torch.tensor

fermionic(A)[source]

Parameters: A (torch.tensor) – an auxiliary field configuration
Returns: \(-\log \det \mathbb{d}\)
Return type: torch.tensor

__call__(A)[source]

Parameters: A (torch.tensor) – an auxiliary field configuration
Returns: \(S(A) =\texttt{S.gauge}(A) + \texttt{S.fermionic}(A) + \texttt{S.normalizing_offset}\).
Return type: torch.tensor

quenched_sample(sample_shape=torch.Size([]))[source]

Provides sample auxiliary fields drawn from the gaussian

\[p(A) \propto \exp\left( - \frac{1}{2} \sum_t A_t (-\Delta t V)^{-1} A_t \right)\]

Parameters: sample_shape (torch.Size) – See the torch.distributions interface for details.
Returns: Shape is [*sample_shape, *spacetime.dims]. Called with no arguments the default is to produce just one sample.
Return type: torch.tensor

set_tuning(ere)[source]

projected(n, s)[source]

Provides a convenience constructor for actions derived from this one that are needed in the implementation of Sector. Shifts the chemical potential and external field, via

\[\begin{align} \mu &\rightarrow \mu + \frac{2\pi i}{2V+1} \frac{n}{\beta} & \vec{h} &\rightarrow \vec{h} + \frac{2\pi i}{2V+1} \frac{2s}{\beta} \hat{h} \end{align}\]

Parameters

n (integer) – Fourier sector for particle number projection
s (half-integer) – Fourier sector for spin projection

Return type

Action

The Grand Canonical Ensemble

With the auxiliary field Action in our hands we can sample the distribution \(\exp(-S)\). An assignment for the auxiliar field is a configuration, and set of configurations is called an ensemble.

We can generate() an ensemble using a sampling scheme (described in Hybrid Monte Carlo) or build one from_configurations() we already have.

tdg.ensemble.GrandCanonical

class tdg.ensemble.GrandCanonical(Action)[source]

Bases: H5able

A grand-canonical ensemble of configurations and associated observables, importance-sampled according to Action.

Note

GrandCanonical also supports a large number of observables.

Parameters: Action (tdg.Action) – An action which describes a Euclidean path integral equal to a Trotterization of the physics of interest.

Action: The action with which the ensemble was constructed.

from_configurations(configurations, weights=None, index=None)[source]

Parameters

configurations (torch.tensor) – A set of pre-computed configurations.
weights (torch.tensor) – Weights for each configuration.
index (torch.tensor) – Where in Markov chain time did each configuration come from.

Returns

the ensemble itself, so that one can do ensemble = GrandCanonical(action).from_configurations(phi).
If weights is None, the weights are all 1. If index is None, the index counts up from 0.

generate(steps, generator, start='hot', progress=<function _no_op>, starting_index=0)[source]

Parameters

steps (int) – Number of configurations to generate.
generator – Something which produces a new configuration if called as generator.step(previous_configuration). Often an hmc instance. May be provided with a default in the future.
start ('hot', 'cold', or torch.tensor) – A hot start begins with a configuration drawn from the quenched action. A cold start beins with the zero configuration. If a tensor is passed that tensor is used as the first configuration.
progress (something which wraps an iterator and provides a progress bar.) – In a script you might use tqdm.tqdm, and in a notebook tqdm.notebook. Defaults to no progress reporting.
starting_index (int) – The generated steps will have [starting_index, starting_index+steps).

Return type

the ensemble itself, so that one can do ensemble = GrandCanonical(action).generate(...).

Populates the index attribute, a torch tensor counting up from starting_index, once for each call to the generator, so that each configuration has an index. This index is kept track of through cut(), every(), binned() (by the Binning).

classmethod continue_from(ensemble, steps, progress=<function _no_op>)[source]

Use the last configuration and generator of ensemble to produce a new ensemble of steps configurations.

Parameters

ensemble (tdg.ensemble.GrandCanonical or an h5py.Group that encodes such an ensemble) – The ensemble to continue. Raises a ValueError if it is not a tdg.ensemble.GrandCanonical or an h5py.Group with an action, generator, and at least one configuration.
steps (int) – Number of configurations to generate.

Returns

An ensemble with steps new configurataions generted in the same way as ensemble.

Return type

tdg.ensemble.GrandCanonical

Todo

The starting weight should automatically be read in; currently not.

measure(*observables)[source]

Compute each @observable in observables; log an error for any unregistered observable.

Parameters: observables (strings) – Names of observables.
Return type: GrandCanonical; itself, now with some observables measured.

Note

If no observables are passed, evaluates every registered @observable.

to_h5(group, _top=True)[source]: Just like tdg.h5.H5able.to_h5() but some data are written as resizable datasets

extend_h5(group, fields=None)[source]

Append the measured values for this ensemble onto the resizable datasets (configurations, index, weights, and observables).

Warning

Currently does not check that the actions match! So sloppy behavior on your part can result in really weird, probably-meaningless data.

Parameters

group (h5 group) – Should be an h5 group written from an ensemble with the same action.
fields – Tuple of properties to extend. If None, extends configurations, index, weights, and all observbles evaluated; does not trigger the measurement of observables.

classmethod from_h5(group, strict=True, observables=None, selection=(), _top=True)[source]

Parameters

group (h5 group) – Where the data for this ensemble is found; same as tdg.h5.H5Data.from_h5().
strict (boole) – Same as tdg.h5.H5Data.from_h5().
observables (iterable of strings) – Which observables should be read. If None read all the observables on disk.
selection (fancy indexing) – A subset of numpy fancy indexing is supported; this selection is used for selecting configurations bassed on their location in the dataset, rather than their .index. Some valid choics are [1,2,3], slice(1,4), slice(1,10,2). If you want a singleton, you should use [1,] or the configurations will have the wrong number of dimensions.

cut(start)[source]

Parameters: start (int) – The number of configurations to drop from the beginning of the ensemble.
Return type: GrandCanonical without some configurations at the start. Useful for performing a thermalization cut.

every(frequency)[source]

Parameters: frequency (int) – The frequency with which to keep configurations.
Return type: GrandCanonical with configurations reduced in size by a factor of the frequency. Useful for Markov Chain decorrelation.

binned(width=1)[source]

Parameters: width (int) – The width of the bins over which to average.
Return type: Binning of the ensemble, with the width specified.

bootstrapped(draws=100)[source]

Parameters: draws (int) – Resamples for uncertainty estimation; see Bootstrap for details.
Return type: Bootstrap built from the ensemble, with the draws specified.