Fix setting of sso.m_ops in heterodyne smesolver and passing through of sc_ops to photocurrent solver.

[hodgestar]

Description
Fix setting of sso.m_ops in heterodyne smesolver and passing through of sc_ops to photocurrent solver.

Related issues or PRs

• Variable used before declaration in smesolve for heterodyne detection. #2078

[hodgestar]

@theodotk Thank you for the bug report. I cherry-picked your fix into this PR. Any thoughts on a simple test case we could add? Ideally we could just pass in custom m_ops in something like test_smesolve_heterodyne but I don't know what the custom m_ops should be, how to calculate the expected result or how to be sure the right m_ops were used.

[theodotk]

@hodgestar Thanks for doing it!
The custom m_ops, as far as I understand, would be of a form [np.sqrt(coef)*(a+a.dag()), -1i*np.sqrt(coef)*(a - a.dag())] where a would be an annihilation operator for a photon in a mode that is detected, and coef represents whatever happens between the emission and the detection (cavity-guide coupling, detection efficiency etc). In the current implementation it would be subset of default m_ops, that are formed as

for c in sso.sc_ops:
    m_ops += [c + c.dag(), -1j * (c - c.dag())]

As for the tests, I can think of something like this

def test_heterodyne_mesolve():
    b = 1  # drive amplitude
    gamma = 1  # spont.  emission rate
    eta = 0.3  # coupling efficiency
    n_steps = 1000
    n_traj = 50
    
    H = np.sqrt(eta*gamma) * b * sigmay()
    c_ops = [np.sqrt(gamma)*sigmam()]
    psi0 = basis(2)
    times = np.linspace(0, np.pi*2, n_steps)

    n_tr = 3

    sme_het = smesolve(
        H,
        psi0,
        times,
        [],
        c_ops,
        e_ops=[sigmax(), sigmay(), sigmaz()],
        store_measurement = True,
        dW_factors=[1e-5, 1e-5],  # to make tests simpler
         method="heterodyne",
         m_ops=[np.sqrt(eta)*sigmax(), np.sqrt(eta)*sigmay()],
         ntraj=n_traj,
        noise=123,  # random seed
    )
    
    assert np.array(sme_het.measurement).shape == (n_traj, n_steps, 1, 2)
    assert all(np.isclose(np.array(sme_het.measurement).mean(axis=0)[:,0,0].T, np.sqrt(eta)*sme_het.expect[0], atol=2e-2))
    
def test_incorrect_m_ops_heterodyne_mesolve():
    
    with pytest.raises(ValueError, match="The measured operators for the heterodyne method supposed to be pairs of quadratures: m_ops should have even length."):
        sme_het_bad_mops = smesolve(
            sigmax(),
            basis(2),
            np.linspace(0,1,10),
            [],
            [sigmam()],
            e_ops=[],
            store_measurement = True,
             method="heterodyne",
             m_ops=[np.sqrt(eta)*sigmax(), np.sqrt(eta)*sigmay(), np.sqrt(eta)*sigmaz()],  # three operators
             ntraj=10,
        )

[hodgestar]

@theodotk Apologies for the delay and thank you for the tests. I've incorporate them now and they look good. It would be nice to make the expected results match a little better. Supplying only atol is problematic when the values themselves are small and currently the relative differences are larger than 1. Any suggestion for making these better without having the test take too long?

[coveralls]

Coverage: 69.957%. Remained the same when pulling 9091a39 on hodgestar:fix/smesolve-4-7-bugfix into 1ce29f7 on qutip:qutip-4.7.X.

[Ericgig]

There should be 2 m_ops per sc_ops with heterodyne. The measurement use the noise of the corresponding sc_ops in it's computation. Having m_ops than effective sc_ops will cause junk measurement.

[theodotk]

Thanks for the PR @hodgestar !
I was just about to write that large rtol values are okay, as at some points of the simulation that values that we are comparing are close to 0, and any noise would be large. I thought that the way to make the atol values smaller would be decreasing dW_factors argument, but it didn't quite work for me.