Calculating two-time correlations for a chain of sites

[BrendonLovett1]

I'd like to be able to calculate two-time correlations in for a chain of sites, each coupled to an environment. This would enable the calculation of effective temperature using the fluctuation-dissipation theorem, as in https://doi.org/10.1103/PhysRevResearch.5.033078.

[gefux]

Hi @BrendonLovett1!
While for single systems there is a compute_correlations() and a compute_correlations_nt() function, for chains of systems this needs to be done by inserting the corresponding control operations (i.e. the left or right hand application of an operator) at the time and cite of interest.

For that you can pass a ChainControl object to the constructor of the PtTebd class. In the ChainControl object you can specify which super operator should be applied at which time step on which site.

I'll follow up with some simple example in the next days to clarify.

[BrendonLovett1]

Hi @gefux Thank you for your help! OK, this sounds very do-able! I'll get back in touch if we make any progress here, and look forward to looking at your example.

[gefux]

We calculate two different two-time correlations, namly the time-ordered $\langle B(t_1+\tau) A(t_1) \rangle$, and the anti time-ordered correlation $\langle B(t_1) A(t_1+\tau) \rangle$, where $A$ and $B$ are operators acting on siteA and siteB respectively.
The time-ordered correlation is computed by evolving the system first to time $t_1$, acting with $A$ to the left of the density matrix, continuing the evolution and then compute the (complex valued) expectation of $B$. The anti time-ordered correlation is computed by also evolving the system first to time $t_1$, but then acting with $B$ to the right of the density matrix, continuing the evolution and then compute the (complex valued) expectation of $A$. As both require the evolution to time $t_1$ we can perform this computation only once and then continue with two separate computations by copying the PtTebd object (or alternatively the augmented mps). We do this in the following.

from copy import deepcopy
import oqupy
import oqupy.operators as op
import numpy as np
import matplotlib.pyplot as plt

# First we set-up a PT-TEBD calculation as is done in the PT-TEBD tutorial ...

N = 5
epsilon = 1.0
J = [1.3, 0.7, 1.2]

system_chain = oqupy.SystemChain(hilbert_space_dimensions=[2]*N)

for n in range(N):
    system_chain.add_site_hamiltonian(
        site=n,
        hamiltonian=0.5*epsilon*op.sigma("z"))

for n in range(N-1):
    for J_gamma, s_gamma in zip(J, [0.5*op.sigma("x"), 0.5*op.sigma("y"), 0.5*op.sigma("z")]):
        system_chain.add_nn_hamiltonian(
            site=n,
            hamiltonian_l=J_gamma*s_gamma,
            hamiltonian_r=s_gamma)

up_dm = op.spin_dm("z+")
down_dm = op.spin_dm("z-")
initial_augmented_mps = oqupy.AugmentedMPS([up_dm, down_dm, down_dm, down_dm, down_dm])

pt_tebd_params = oqupy.PtTebdParameters(
    dt=0.2,
    order=2,
    epsrel=1.0e-6)

pt_tebd = oqupy.PtTebd(
    initial_augmented_mps=initial_augmented_mps,
    system_chain=system_chain,
    process_tensors=[None, None, None, None, None],
    parameters=pt_tebd_params,
    dynamics_sites=[0, 1, 2, 3, 4])


# Now let's define the two time correlations we'd like to calculate:

num_steps = 20
t1 = 10              # the 10th time step 
opA = op.sigma("x")
opB = op.sigma("x")
siteA = 0
siteB = 1


# We first compute the evolution until time just before $t_1$:

pt_tebd.compute(t1-1, progress_type="bar");


# Then we create the control objects and enter them to two copies of the pt_tebd object:

chain_control_leftA = oqupy.ChainControl(hilbert_space_dimensions=[2]*N)
chain_control_leftA.add_single_site_control(control=op.left_super(opA), site=siteA, step=t1)
chain_control_rightB = oqupy.ChainControl(hilbert_space_dimensions=[2]*N)
chain_control_rightB.add_single_site_control(control=op.right_super(opB), site=siteB, step=t1)

pt_tebd_leftA = pt_tebd
pt_tebd_rightB = deepcopy(pt_tebd)
pt_tebd_leftA.chain_control = chain_control_leftA
pt_tebd_rightB.chain_control = chain_control_rightB


# We can then compute the dynamics for both cases, and extract the correlation as a function of the second time.

res_leftA = pt_tebd_leftA.compute(num_steps)
res_rightB = pt_tebd_rightB.compute(num_steps)

ts, corr_time_ordered = res_leftA['dynamics'][siteB].expectations(opB, real=False)
ts, corr_anti_ordered = res_rightB['dynamics'][siteA].expectations(opA, real=False)
tau = ts[t1:] - ts[t1]


# Finally here we plot real and imaginary parts of the correlations:

plt.plot(tau, corr_time_ordered.real[t1:],
         color="C0", linestyle="solid",
         label=r"$\langle B(t_1+\tau) A(t_1) \rangle$")
plt.plot(tau, corr_time_ordered.imag[t1:],
         color="C0", linestyle="dotted")
plt.plot(tau, corr_anti_ordered.real[t1:],
         color="C1", linestyle="solid",
         label=r"$\langle B(t_1) A(t_1+\tau) \rangle$")
plt.plot(tau, corr_anti_ordered.imag[t1:],
         color="C1", linestyle="dotted")
plt.xlabel(r"$\tau$")
plt.ylabel("correlation")
plt.legend()

plt.show()

Let us know if this answers your question.