In the previous tutorial, we showed how to use `pymf` to solve a simple 1D Hubbard model with onsite interactions.
In this tutorial, we will apply `pymf` to more complex system: graphene with onsite $U$ and nearest-neighbour $V$ interactions.
The system is more complicated in every aspect: the lattice structure, dimension of the problem, complexity of the interactions.
And yet, the workflow is the same as in the previous tutorial and remains simple and straightforward.
## Building the system with `kwant`
### Non-interacting part
As in the previous tutorial, we could construct a tight-binding dictionary of graphene by hand, but instead it is much easier to use [`kwant`](https://kwant-project.org/) to build the system.
For a more detailed explanation on `kwant` see the [tutorial](https://kwant-project.org/doc/1/tutorial/graphene).
The `bulk_graphene` object is a `kwant.Builder` object that represents the non-interacting graphene system.
To convert it to a tight-binding dictionary, we use the {autolink}`~pymf.kwant_helper.utils.builder_to_tb` function:
```{code-cell} ipython3
h_0 = utils.builder_to_tb(bulk_graphene)
```
### Interacting part
We utilize `kwant` to build the interaction tight-binding dictionary as well.
To define the interactions, we need to specify two functions:
*`onsite_int(site)`: returns the onsite interaction matrix.
*`nn_int(site1, site2)`: returns the interaction matrix between `site1` and `site2`.
We feed these functions to the {autolink}`~pymf.kwant_helper.utils.build_interacting_syst` function, which constructs the `kwant.Builder` object encoding the interactions.
All we need to do is to convert this object to a tight-binding dictionary using the {autolink}`~pymf.kwant_helper.utils.builder_to_tb` function.
```{code-cell} ipython3
def onsite_int(site, U):
return U * sx
def nn_int(site1, site2, V):
return V * np.ones((2, 2))
builder_int = utils.build_interacting_syst(
builder=bulk_graphene,
lattice=graphene,
func_onsite=onsite_int,
func_hop=nn_int,
max_neighbor=1
)
params = dict(U=0.2, V=1.2)
h_int = utils.builder_to_tb(builder_int, params)
```
Because `nn_int` function returns the same interaction matrix for all site pairs, we set `max_neighbor=1` to ensure that the interaction only extends to nearest-neighbours and is zero for longer distances.
## Computing expectation values
As before, we construct {autolink}`~pymf.model.Model` object to represent the full system to be solved via the mean-field approximation.
We then generate a random guess for the mean-field solution and solve the system:
```{code-cell} ipython3
filling = 2
model = pymf.Model(h_0, h_int, filling=2)
int_keys = frozenset(h_int)
ndof = len(list(h_0.values())[0])
guess = pymf.generate_guess(int_keys, ndof)
mf_sol = pymf.solver(model, guess, nk=18)
h_full = pymf.add_tb(h_0, mf_sol)
```
To investigate the effects of interaction on systems with more than one degree of freedom, it is more useful to consider the expectation values of various operators which serve as order parameters.
For example, we can compute the charge density wave (CDW) order parameter which is defined as the difference in the charge density between the two sublattices.
To calculate operator expectation values, we first need to construct the density matrix via the {autolink}`~pymf.mf.construct_density_matrix` function.
We then feed it into {autolink}`~pymf.observables.expectation_value` function together with the operator we want to measure.
In this case, we compute the CDW order parameter by measuring the expectation value of the $\sigma_z$ operator acting on the graphene sublattice degree of freedom.
This phase diagram has gap openings at the same places as shown in the [literature](https://arxiv.org/abs/1204.4531).
We can now use the stored results in `mf_sols` to fully map out the phase diagram with order parameters.
On top of the charge density wave (CDW), we also expect a spin density wave (SDW) in a different region of the phase diagram.
We construct the SDW order parameter with the same steps as before, but now we need to sum over the expectation values of the three Pauli matrices to account for the $SU(2)$ spin-rotation symmetry.
Kostas Vilkelis authored