18 fix dissipator coupling bug

[EthanMcGarrigle-SNL]

Fixed a bug where coupling operators would only store an upper half or lower half of the matrix because of assumed usage in a Hamiltonian class. This led to Lindblad operators being created incorrectly if the Lindblad operator had content in both the upper and lower portions of its matrix.

Changes:

1. Coupling operators parse and store all couplings.

2. Coupling and General operators have an attribute "hermitian_couplings" which will return the unique couplings. For a General Operator with no couplings, None is returned.

3. The Hamiltonian class manages coupling operators by working from the hermitian couplings and all previous logic executes as before. i.e. a Coupling operator's Hamiltonian contribution is built via its static matrix * rate + the hermitian conjugate.

Note: During testing of different Lindblad operator choices, if the user specifies a \sigma_X type Lindblad operator in the lab frame that is highly off-resonant with the qubit, then required robust numerical integration methods are required since the transformed Lindblad operator is highly-oscillatory in the interaction picture. Therefore, there are several changes I made to use the 'zvode' integration method:

4. Custom math "zvode" method required some syntax fixes for attributes.

5. Made "atol" and "rtol" zvode class attributes so they could set by the user upstream as keyword arguments. The default atol and rtol values are too aggressive (near machine precision!) and lead to excessively slow integration behavior for problems that don't need such tight tolerances.

6. Added some parsing to State to return both the times vector and the states list from "propagate_using_master_equation()" when using 'zvode' since 'zvode' uses an adaptive time stepper and alters the list of times actually used in the simulation. i.e. If a user specifies a "times" vector, the state list from 'zvode' does not match the user's times list, so 'zvode' should give the user the correct times vector determined during numerical integration. EDIT 07/27/2026, instead of this approach, we have fixed 'zvode' to always return the solution at the exact time points as requested by the user. This is better than having heterogeneous output from the "propagate_using...()" method in state.

7. Added a test to test_dissipators.py to test the specific example of a qubit with an interaction-frame Lindblad operator L = \sqrt{gamma} \sigma_X. It tests both integration and agreement of the final |0> population with a reference.

[brandonruzic]

Re: "Coupling and General operators have an attribute "hermitian_couplings" which will return the unique couplings," why wouldn't we just fix the Couplings to store all off-diagonal elements? Isn't it confusing to have two separate attributes?

[brandonruzic]

Re: "Added some parsing to State to return both the times vector and the states list from "propagate_using_master_equation()" when using 'zvode' since 'zvode' uses an adaptive time stepper and alters the list of times actually used in the simulation. i.e. If a user specifies a "times" vector, the state list from 'zvode' does not match the user's times list, so 'zvode' should give the user the correct times vector determined during numerical integration,"

I believe that propagate_using_master_equation already did return both the time array and the density matrices on that time grid.-- i.e., (times, rhos) -- and these times are not the ones used for numerical integration. The time grid specified, "time-evals" is meant to be the time grid on which the results are returned, independent of the integration time grid. Was this not the functionality of the zvode integration method?

[EthanMcGarrigle-SNL]

Re: "Coupling and General operators have an attribute "hermitian_couplings" which will return the unique couplings," why wouldn't we just fix the Couplings to store all off-diagonal elements? Isn't it confusing to have two separate attributes?

The Coupling Operators and General Operators do store all off-diagonal elements as a property "couplings", which returns a list of coupling elements.

Having a separate attribute called "hermitian_couplings" helped give the Hamiltonian class easy access to the unique couplings it uses to generate the "Hint" contributions in the current implementation, but this can be done without needing "hermitian_couplings" to be its own property attribute. I'll remove the "hermitian_couplings" attribute since it's essentially a wrapper anyway and doesn't get used anywhere else besides the Hamiltonian; this should be reflected in the latest commit.

[EthanMcGarrigle-SNL]

I believe that propagate_using_master_equation already did return both the time array and the density matrices on that time grid.-- i.e., (times, rhos) -- and these times are not the ones used for numerical integration. The time grid specified, "time-evals" is meant to be the time grid on which the results are returned, independent of the integration time grid. Was this not the functionality of the zvode integration method?

The "propagate_using_master_equation" method does not return any time information. For reference, see a screenshot below from "propagate_using_master_equation()" from the main branch:

I agree that zvode should return the states corresponding to the "time_evals" specified by the user, but I think the answer is "Yes", i.e. zvode can return a different time grid than the one specified by the user, but they are very close. Depending on how the user specifies their time grid, zvode's time grid may differ by an extra element, so that if a user gives 20 time points, zvode's time grid has 21 timepoints. This doesn't happen all the time.

An example:

1. times = np.linspace(0, duration, int(duration/dt) + 1) will correspond exactly with the zvode times output, so zvode's list of States would have the same length as "times".

2. times = np.linspace(0., duration, N_points) may be off by 1 point with the times used in zvode, so the list of States will be 1 element longer than "times". This happens only for some choices of "N_points", so it's inconsistent.

I think the user should be able to specify the time grid using either way.

Please see the latest commit: I've made a change to zvode's method in custom_math.py to ensure the states are evaluated at the corresponding time evals points exactly as specified by the user. I've updated the "propagate_using" method accordingly, so it doesn't return different objects if using zvode. In the zvode custom_math method, I've left the previous code commented out in case you want to revert back to the old method, which used a while loop instead of an explicit for loop over the "time_evals". If you think my change is okay, I will remove the commented code.