rydiqule.timesolvers.solve_time

rydiqule.timesolvers.solve_time(sensor: Sensor, end_time: float, num_pts: int, init_cond: ndarray | None = None, doppler: bool = False, doppler_mesh_method: UniformMethod | IsoPopMethod | SplitMethod | DirectMethod | None = None, sum_doppler: bool = True, weight_doppler: bool = True, n_slices: int | None = None, solver: Callable | Literal['scipy', 'cyrk'] = 'scipy', **kwargs) → Solution

Solves the response of the optical sensor in the time domain given the its time-dependent inputs

If insuffucent system memory is available to solve the system all at once, system is broken into “slices” of manageable memory footprint which are solved indivudually. This slicing behavior does not affect the result. All couplings that include a “time_dependence” argument will be solved in the time domain.

A number of solver backends work with rydiqule, but the default "scipy" ivp solver is the is recommended backend in almost all cases, as it is the most fully-featured and documented. Advanced users have the ablity to define their own solver backends by creating a function that follows the call signature for rydiqule timesolver backends. Additional arguments to the solver backend can be supplied with **kwargs.

Parameters:

• sensor (Sensor) – The sensor object representing the atomic/laser arrangement of the system.

• end_time (float) – Amount of time, in microseconds, for which to simulate the system

• num_pts (int) – The number of points along the range (0, end_time) for which the solution is evaluated. This does not affect the number of funtion evaluations during the solve, rather the spacing of the points in the reported solution.

• init_cond (numpy.ndarray or None, optional) – Density matrix representing the initial state of the system. If specified, the shape should be either (n) in the case of a single initial condition for all parameter values, or should be of shape (*l, n) matching the output shape of a steady state solve if the initial condition may be different for different combinations of parameters. If None, will solve the problem in the steady state with all time-dependent fields at their \(t=0\) value and use the solution as the initial condition. Other possible manual options might include a matrix populated by zeros representing the entire population in the ground state. Defaults to None.

• doppler (bool, optional) – Whether to account for doppler shift among moving atoms in the gas. If True, the solver will implicitly define a velocity distribution for particles in the cell, solve the problem for each velocity class, and return a weighted average of the results. Note that solving in this manner carries a substantial performance penalty, as each doppler velocity class is solved as its own problem. If solved with doppler, only axis specified by a "kvec"' argument in one of the sensor couplings will be average over. The time solver currently supports doppler averaging in any number of spatial dimensions, up to the limit of 3 imposed by the macroscopic physical world. Defaults to `False.

• doppler_mesh_method (dict, optional) – Dictionary that controls the doppler meshing method. Exact details of this are found in the documentation of doppler_classes(). Ignored if doppler=False. Default is None.

• sum_doppler (bool, optional) – Whether to average over doppler classes after the solve is complete. Setting to false will not perform the sum, allowing viewing of the weighted results of the solve for each doppler class. Ignored if doppler=False. Default is True.

• weight_doppler (bool) – Whether to apply weights to doppler solution to perform averaging. If False, will not apply weights or perform a doppler_average, regardless of the value of sum_doppler. Changing from default intended only for internal use. Ignored if doppler=False or sum_doppler=False. Default is True.

• n_slices (int or None, optional) – How many sets of equations to break the full equations into. The actual number of slices will be the largest between this value and the minumum number of slices to solve the system without a memory error. If None, solver uses the minimum number of slices required to solve without a memoryError. Defaults to None.

• solver ({"scipy", "cyrk"} or callable) –

  The backend solver used to solve the ivp generated by the sensor. All string values correspond to backend solvers built in to rydiqule. Valid string values are:

  • ”scipy”: Solves equations with scipy.integrate.solve_ivp(). The default, most stable, and well-supported option.

  • ”cyrk”: Solves jit-compiled equations with a cython compiled RK solver from CyRK. Due to some jit compilation, only faster for moderate length problems (ie problems with a moderate number of required time steps).

  Additionally, can be specified with a callable that matches rydiqule’s time-solver convention, enabling using a custom solver backend.

  Note

  Unless otherwise noted, backends other than scipy are considered experimental. Issues with their use are considered features not fully implemented rather than bugs.

• **kwargs (Additional keyword arguments passed to the backend solver.) – See documentation of the relevant solver (i.e. scipy.integrate.solve_ivp()) for details and supported arguments.

Returns:

An object contining the solution and related information. Timesolver-specific defined attributes are t and init_cond, corresponding respectively to the times at which the solution is sampled and the initial conditions used for the solve.

Return type:

Solution

Examples

A basic solve for a 3-level system would have a “density matrix” solution of size 8 (3^2-1). Here we use a trivial time dependence for demonstration purposes, but in practice the time dependence is likely more complicated. Below the most basic use of solve_time is demonstrated

>>> s = rq.Sensor(3)
>>> td = lambda t: 1
>>> s.add_coupling((0,1), detuning = 1, rabi_frequency=1)
>>> s.add_coupling((1,2), detuning = 2, rabi_frequency=2, time_dependence=td)
>>> s.add_transit_broadening(0.1)
>>> end_time = 10 #microseconds
>>> n_pts = 1000 #interpoleted points in solution
>>> sol = rq.solve_time(s, end_time, n_pts)
>>> print(type(sol))
<class 'rydiqule.sensor_solution.Solution'>
>>> print(type(sol.rho))
<class 'numpy.ndarray'>
>>> print(sol.rho.shape)
(1000, 8)

Defining an array-like parameter will automatically calculate the density matrix solution for every value. Here we use 11 values, resulting in 11 density matrices. The axis_labels attribute of the solution can clarify which axes are which.

>>> s = rq.Sensor(3)
>>> td = lambda t: 1
>>> det = np.linspace(-1,1,11)
>>> s.add_coupling((0,1), detuning = det, rabi_frequency=1)
>>> s.add_coupling((1,2), detuning = 2, rabi_frequency=2, time_dependence=td)
>>> s.add_transit_broadening(0.1)
>>> sol = rq.solve_time(s, end_time, n_pts)
>>> print(sol.rho.shape)
(11, 1000, 8)
>>> print(sol.axis_labels)
['(0,1)_detuning', 'time', 'density_matrix']

As expected, multiple axes of scanned parameters are handled the same way as they are in the steady-state case, with the expected additions from the time solver.

>>> s = rq.Sensor(3)
>>> td = lambda t: 1
>>> det = np.linspace(-1,1,11)
>>> s.add_coupling((0,1), detuning = det, rabi_frequency=1)
>>> s.add_coupling((1,2), detuning = det, rabi_frequency=2, time_dependence=td)
>>> s.add_transit_broadening(0.1)
>>> sol = rq.solve_time(s, end_time, n_pts)
>>> print(sol.rho.shape)
(11, 11, 1000, 8)
>>> print(sol.axis_labels)
['(0,1)_detuning', '(1,2)_detuning', 'time', 'density_matrix']

If the solve uses doppler broadening, all doppler classes will be computed a weighted average will be taken over the doppler axis, and the shape of the solution will not change. While this is the desired behavior in most situations, the sum_doppler argument can be used to override this behavior, leave (or more) solution axis corresponding to different doppler classes.

>>> s = rq.Sensor(3, vP=1)
>>> td = lambda t: 1
>>> det = np.linspace(-1,1,11)
>>> s.add_coupling((0,1), detuning = det, rabi_frequency=1, kvec=(4,0,0))
>>> s.add_coupling((1,2), detuning = 2, rabi_frequency=2, kvec=(-4,0,0), time_dependence=td)
>>> s.add_transit_broadening(0.1)
>>> end_time = 10 #microseconds
>>> n_pts = 1000 #interpoleted points in solution
>>> sol_dop = rq.solve_time(s, end_time, n_pts, doppler=True)
>>> sol_dop_nosum = rq.solve_time(s, end_time, n_pts, doppler=True, sum_doppler=False)
>>> print(sol_dop.rho.shape)
(11, 1000, 8)
>>> print(sol_dop_nosum.rho.shape)
(561, 11, 1000, 8)
>>> print(sol_dop.axis_labels)
['(0,1)_detuning', 'time', 'density_matrix']
>>> print(sol_dop_nosum.axis_labels)
['doppler_0', '(0,1)_detuning', 'time', 'density_matrix']