Time-Dependent Rabi Frequency and Detuning

[ComputationalBasis]

Hello everyone,

I am currently trying to obtain a QPT phase diagram of a chain of 8 atoms. This articles does it https://doi.org/10.1038/s41586-024-08353-5 (fig. 1.b) using a square lattice of 16x16 atoms. I noticed that to observe the QPT, one has to make the Rabi frequency and the detuning time dependent (from negative values to positive values). In this article, https://doi.org/10.1038/s41586-019-1070-1 (fig 1.c. inset) they do exactly what I am trying to recreate, but with a chain of ~50 atoms.

I do not have any problem with time-independent Hamiltonians, but I am having trouble to simulate the dynamics when Omega and delta are time-dependent. I tried using Heaviside functions to define a similar time-dependent that is being used in those articles, but I was always getting an error. I have also read that the phase boundaries are calculated using Density Matrix Renomarlization Group techniques, and I also do not know if I also have to use it in order to get a phase diagram, as it is something completely new for me.

The motivation came from this Thesis I found (page 77) Lucas Leclerc. Quantum computing with Rydberg atoms: control and modelling for quantum simulation and practical algorithms. Optics [physics.optics]. Université Paris-Saclay, 2024. English. �NNT:
2024UPASP046�. �tel-04745992�, but I am not really sure if the QPT is obtainable using Qadence, or maybe it is better to study another topic (maybe something simpler), if anyone has a better option / approach using Rydberg atoms, please let me know.

I would paste my code here but I do not think it is worth it, I am more interested in how to obtain the desired time-dependent on Omega and delta (if it is even possible) and obtain the wave function when evolving the system under this time-dependent Hamiltonian.

I was told that @vytautas-a could help me with time-dependent questions.

[jpmoutinho]

Hey @ComputationalBasis!

Could you post a small script of what you are running and the error you are getting?

Note that analytical Heaviside functions are not differentiable, so that could be one cause for errors.

[ComputationalBasis]

Thanks, Joao.

The logistic function is exactly what I was using, at least for a first approximation to my desired result. I do not think I need such control over the pulses, I was just wondering if there was a way of doing it simple. Also, I was not aware of that Pulser documentation, that seems very useful. I will keep trying things and see if I can finally get the QPT diagram.

[jpmoutinho]

Pulser is a separate package all-together, it implements the Rydberg physics directly. Qadence has the features for defining more general Hamiltonians, but it is not as polished (yet) when it comes to defining different pulse shapes. What you are doing is definitely interesting to the Qadence development. In theory, approximating square pulses with combinations of logistic functions should be possible, but we have not tested it in detail. If you come up with nice implementations, we would be very happy to have you contribute them to the documentation, if you are interested.

[ComputationalBasis]

So it is a separate package, got it. Still, I want to stick with Qadence, as I want to use a Digital-Analog approach with Rydberg atoms. I will then try to approximate the pulses in terms of a Fourier expansion, and see if I get the same (or at least something approximate) result of the adiabatic preparation of an antiferromagnetic state. If I manage to get something promising, I will for sure let you know. Having my little contributtion to the Qadence documentation would be really cool.

Thanks for everything, Joao, you really are a good dev :)

[RolandMacDoland]

Hi @ComputationalBasis if you find any of the above answers useful, could you please mark them ? That'd be helpful for us to monitor which questions have been addressed :) Thanks.
PS. Same would hold for the comment below.

[ComputationalBasis]

Yes, of course. I can only mark one, but thank you all for the help, you are all very helpful :). The code finally computed everything with those suggestions but I need to work a little on something. If I have more questions I will let you all know 👍

[ComputationalBasis]

Hello,

I wrote the code for the adiabatic preparation of the antiferromagnetic state with n atoms, but I am always getting a divergence errror, it says the following:

"Maximum number of time steps reached in adaptive time step ODE solver at time t = 0.015 ('max_steps=10000'). This is likely due to a diverging solution. Try increasing the maximum number of steps, or use a different solver.

I already tried to increase the steps, but it always gives the error at the same time t. I also tried decreasing the number of atoms in the register from 8 to 4, but that didn't solve the problem. I did not try using a different solver, if anyone knows what solver should I use here, I would really appreciate it (maybe it fixes the problem). The code I used is:

import torch
from torch import tensor
from qadence import TimeParameter, FeatureParameter, X, Y, Z, N, add, kron, Register, HamEvo, run, BackendName, PI, expectation
from qadence.analog.constants import C6_DICT
from pyqtorch.utils import SolverType
import numpy as np
from sympy import cos, sin, exp

ode_solver = SolverType.DP5_SE  # time-dependent Schrodinger equation solver method
n_steps_hevo = 10000  # integration time steps used by solver

# Parameters
c_6 = C6_DICT[60] # [rad*μm^6/μs]
U = 2*PI # [rad/μs]
omega_max = 2*U
delta_max = 2*U

# Register
n = 6
dx = (c_6 / U) ** (1/6)  # μm
reg = Register.circle(n_qubits=n, spacing=dx)


# Initial state: ground state
initial_state = torch.zeros(2**n, dtype=torch.complex128).unsqueeze(0)
initial_state[0] = 1

# Hamiltonian parameters

t = TimeParameter("t")
omega_param = FeatureParameter("omega") # [rad/μs]
delta_param = FeatureParameter("delta") # [rad/μs]

a = 500
omega_t = 2 * t * (1/(1+exp(-a*t)) - 1/(1+exp(-a*(t-1/2)))) + 1/(1+exp(-a*(t-1/2))) - 1/(1+exp(-a*(t-7/2))) + (8-2*t) * 1/(1+exp(-a*(t-7/2))) - 1/(1+exp(-a*(t-4)));
delta_t = -(1/(1+exp(-a*t)) - 1/(1+exp(-a*(t-1/2)))) + (-4/3 + 2*t/3) * (1/(1+exp(-a*(t-1/2))) - 1/(1+exp(-a*(t-7/2)))) + 1/(1+exp(-a*(t-7/2))) - 1/(1+exp(-a*(t-4)));


Rb = (c_6/omega_param)**(1/6) 
cociente = Rb/dx

# Hamiltonian
h_x = (omega_param * omega_t / 2) * add(X(i) for i in range(n))
h_n = -1.0 * delta_param * delta_t * add(N(i) for i in range(n))

h_d = h_x + h_n  # Driving Hamiltonian

h_nn = U * add(N(i) @ N((i+1)%n) for i in range(n)) # Ising

h = h_d + h_nn

# Duration time evolution
duration = 4 * 1e3  #  [ns] 1e3 ns = 1 μs


# Omega and delta values
param_values = {
    "omega": tensor(omega_max),
    "delta": tensor(delta_max)
}

# Evolution
evo = HamEvo(h, parameter=t, duration=duration/1000)
wf = run(reg, evo, state=initial_state, values=param_values, backend=BackendName.PYQTORCH)

# Sample
xs = sample(reg, evo, state=initial_state, values=param_values, n_shots=10000)

Here you can see the time dependent explicit form of Omega and delta that I used, I think that it covers a fairly accurate representation of the desired form of the pulses. When increasing the parameter a, the pulse form is more and more accurate. The code computational runtime is also pretty big for me, it is like 1 hour. This is also a problem to be fixed which I don't really know how but it is something secundary. Another minor problem is a warning I am getting everytime I run the code, which says:

Warning (from warnings module): File "C:\Users\Usuario\AppData\Local\Programs\Python\Python312\Lib\site-packages\pyqtorch\embed.py", line 133 self.arraylike_fn(symbol_or_numeric, device=self.device) UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).

[jpmoutinho]

Thanks for the script @ComputationalBasis, overall it looks good :)

• Note that buy running run() and then sample() you are actually doing the simulation work twice. If you wish to both save the wavefunction but also sample it, it's best to only to run(), and then implement your own sampling routine from the probabilities.
• The warning you are getting seems like it would be something on our end, but it's not critical.
• The runtime is long, did you try with the other solvers? @vytautas-a knows best that I think the Krylov solver is supposed to be faster.

[vytautas-a]

Thanks for the code @ComputationalBasis. I would suggest experimenting a bit with waveform definition for amplitude and detuning. From my experience creating such "pulse" with explicit exponentials is not very solver-fiendly, meaning it might require more time to converge. Try to use tanh functions instead of exp. Also as @jpmoutinho suggested, try using the KRYLOV_SE solver.

[ComputationalBasis]

Thank you both for the suggestions :) I will try everything and see if now the code improves. I also commented the run line, you are right that the computational time is increased for no reason if I want to sample the final state.