diff --git a/ci_env.yml b/ci_env.yml
index 22d198f..1b17e51 100644
--- a/ci_env.yml
+++ b/ci_env.yml
@@ -4,7 +4,7 @@ channels:
 dependencies:
   - libblas =*=*mkl
   - adcc >=0.16.0
-  - respondo >=0.0.5
+  - respondo >=0.0.6
   - numpy >=1.14
   - isort
   - ruff
diff --git a/examples/cme.py b/examples/cme.py
new file mode 100644
index 0000000..28fa012
--- /dev/null
+++ b/examples/cme.py
@@ -0,0 +1,127 @@
+"""
+Compute the Cotton-Mouton constant according to Eq. 2 in 10.1021/acs.jpca.3c04963
+"""
+import adcc
+import numpy as np
+from pyscf import gto, scf
+from scipy import constants
+
+from responsefun import evaluate_property_isr, TransitionMoment
+from responsefun.symbols_and_labels import (O, n, p, m, mu_a, mu_b, m_a, m_b, m_c, m_d, xi_cd,
+w_n, w_p, w_o, w_2, w_m, w_3, w, w_1
+)
+from responsefun.AdccProperties import DiamagneticMagnetizability
+
+# run SCF in PySCF
+mol = gto.M(
+    atom="""
+    O 0 0 0
+    H 0 0 1.795239827225189
+    H 1.693194615993441 0 -0.599043184453037
+    """,
+    unit="Bohr",
+    basis="sto-3g"
+)
+scfres = scf.RHF(mol)
+scfres.kernel()
+
+# run ADC calculation using adcc
+state = adcc.adc2(scfres, n_singlets=1)
+
+# define (first term of) symbolic SOS expressions
+alpha_sos_expr = (
+        TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, O) / (w_n - w)
+        + TransitionMoment(O, mu_b, n) * TransitionMoment(n, mu_a, O) / (w_n + w)
+)
+xi_para_sos_expr = (
+        TransitionMoment(O, m_a, n) * TransitionMoment(n, m_b, O) / (w_n - w)
+        + TransitionMoment(O, m_b, n) * TransitionMoment(n, m_a, O) / (w_n + w)
+)
+eta_dia_sos_expr = (
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, p, shifted=True)
+    * TransitionMoment(p, xi_cd, O) / ((w_n - w_o) * (w_p - w_2))
+)
+eta_para_term_I_sos_expr = (
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, m, shifted=True)
+    * TransitionMoment(m, m_c, p, shifted=True) * TransitionMoment(p, m_d, O)
+    / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
+)
+eta_para_term_II_sos_expr = (
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, O)
+    * TransitionMoment(O, m_c, m) * TransitionMoment(m, m_d, O)
+    / ((w_n - w_o) * (w_m - w_3) * (w_m + w_2))
+)
+
+w_ruby = 0.072
+gauge_origin = "mass_center"
+
+# calculate tensors
+alpha_tens = evaluate_property_isr(
+    state, # ExcitedStates object returned by adcc calculation
+    alpha_sos_expr, # symbolic SOS expression
+    [n], # indices of summation
+    excluded_states=O,# states excluded from summation (here: ground state)
+    freqs_in=(w, w_ruby), # incoming frequencies
+    freqs_out=(w, w_ruby), # outcoming frequencies
+    conv_tol=1e-4, # convergence tolerance for response solver
+    gauge_origin=gauge_origin # gauge origin for operator integrals
+)
+
+xi_para_tens = evaluate_property_isr(
+    state, xi_para_sos_expr, [n], excluded_states=O,
+    freqs_in=(w, 0), freqs_out=(w, 0),
+    conv_tol=1e-4, 
+    gauge_origin=gauge_origin
+)
+
+xi_dia_tens = DiamagneticMagnetizability(state, gauge_origin=gauge_origin).gs_moment
+
+eta_dia_tens = evaluate_property_isr(
+    state, eta_dia_sos_expr, [n, p],
+    perm_pairs=[(mu_a, -w_o), (mu_b, w_1), (xi_cd, w_2)], # pairs to be permuted
+    freqs_in=[(w_1, w_ruby), (w_2, 0)],
+    freqs_out=(w_o, w_1+w_2),
+    excluded_states=O,
+    conv_tol=1e-4,
+    gauge_origin=gauge_origin
+)
+
+eta_para_tens_I = evaluate_property_isr(
+    state, eta_para_term_I_sos_expr, [n, m, p],
+    perm_pairs=[(mu_a, -w_o), (mu_b, w_1), (m_c, w_2), (m_d, w_3)],
+    excluded_states=O,
+    freqs_in=[(w_1, w_ruby), (w_2, 0), (w_3, 0)],
+    freqs_out=(w_o, w_1+w_2+w_3),
+    conv_tol=1e-4,
+    gauge_origin=gauge_origin
+)
+
+eta_para_tens_II = evaluate_property_isr(
+    state, eta_para_term_II_sos_expr, [n, m],
+    perm_pairs=[(mu_a, -w_o), (mu_b, w_1), (m_c, w_2), (m_d, w_3)],
+    excluded_states=O,
+    freqs_in=[(w_1, w_ruby), (w_2, 0), (w_3, 0)],
+    freqs_out=(w_o, w_1+w_2+w_3),
+    conv_tol=1e-4,
+    gauge_origin=gauge_origin
+)
+
+# add dia- and paramagnetic contributions
+xi_tot = xi_dia_tens + xi_para_tens
+eta_tot = eta_dia_tens + (eta_para_tens_I - eta_para_tens_II)
+
+# calculate anisotropies
+alpha_aniso = alpha_tens[2, 2] - alpha_tens[0, 0]
+xi_aniso = xi_tot[2, 2] - xi_tot[0, 0]
+eta_aniso = 1/15 * (7 * eta_tot[0, 0, 0, 0] - 5 * eta_tot[0, 0, 1, 1] + 2 * eta_tot[2, 2, 2, 2]
+                   - 2 * eta_tot[0, 0, 2, 2] - 2 * eta_tot[2, 2, 0, 0] + 12 * eta_tot[0, 2, 0, 2])
+
+# compute CME constant in cgs system [cm^3 G^–2 mol^–1] at T=298.15K 
+pi = constants.pi
+N_A = constants.N_A
+Hartree = constants.physical_constants['atomic unit of energy'][0]
+k_B = constants.k/Hartree
+T = 273.15
+const_CME = ((2 * np.pi *N_A)/27 * (eta_aniso + 2/ (15 * k_B * T ) * alpha_aniso * xi_aniso))
+const_CME_cgs = const_CME * 2.68211e-44 * 1e20
+print(f"The Cotton-Mouton constant at 298.15 K is {const_CME_cgs:.2f} (10^-20 cm^3 G^–2 mol^–1).")
diff --git a/examples/tpcd.py b/examples/tpcd.py
new file mode 100644
index 0000000..81c0fb0
--- /dev/null
+++ b/examples/tpcd.py
@@ -0,0 +1,101 @@
+"""
+Compute the two-photon circular dichroism rotatory strength in the velocity gauge
+according to 10.1021/acs.jpca.5c02108 eqs. 21-27.
+"""
+import adcc
+import numpy as np
+from pyscf import gto, scf
+
+from responsefun import evaluate_property_isr, TransitionMoment
+from responsefun.misc import epsilon
+from responsefun.symbols_and_labels import (O, f, n, m_b, mup_a, mup_b, mup_c, 
+                                            qp_ab, w_a, w_b, w_f, w_n)
+
+# The calculation of two-photon circular dichroism requires three different two-photon tensors.
+# To avoid code duplication, the following two functions are defined.
+def compute_tp_tensor(state, sos_expr, n_f, perm_pairs, gauge_origin, conv_tol=1e-4):
+    tensor = evaluate_property_isr(
+        state,  # ExcitedStates object returned by adcc calculation
+        sos_expr,  # first term of symbolic SOS expression
+        [n],  # indices of summation
+        excluded_states=O,  # states excluded from summation (here: ground state)
+        excited_state=n_f,  # excited state of interest (here: final state)
+        freqs_in=[(w_a, w_f / 2), (w_b, w_f / 2)],  # incoming frequencies
+        perm_pairs=perm_pairs,  # pairs to be permuted
+        gauge_origin=gauge_origin,  # gauge origin for operator integrals
+        conv_tol=conv_tol  # convergence tolerance for response solver
+    )
+    return tensor
+
+# b1, b2, and b3 are defined by the experimental setup.
+def compute_rotatory_strength(tpcd_data, final_state, b1=6, b2=2, b3=-2):
+    external_energy = tpcd_data["excitation_energy_uncorrected"][final_state]/2
+    Spp = tpcd_data[f"state_{final_state}"]["Spp"]
+    Qpp = tpcd_data[f"state_{final_state}"]["Qpp"]
+    Mp = tpcd_data[f"state_{final_state}"]["Mp"]
+    B1_p = 1/(external_energy**3) * np.einsum("ps,ps->", np.conjugate(Mp), Spp)
+    B2_p = 1/(2 * (external_energy**3)) * np.einsum("ps,ps->", np.conjugate(Qpp), Spp)
+    B3_p = 1/(external_energy**3) * np.einsum("ss->", np.conjugate(Mp)) * np.einsum("pp->", Spp)
+    R_TP_p = 1.0 * b1 * B1_p + b2 * B2_p + 1.0 * b3 * B3_p
+    return R_TP_p
+
+# run SCF in PySCF
+mol = gto.M(
+    atom="""
+    O       0.000000     0.000000     0.000000
+    O       1.480000     0.000000     0.000000
+    H      -0.316648     0.000000     0.895675
+    H       1.796648     0.775678    -0.447838
+    """,
+    unit="Angstrom",
+    basis="sto-3g",
+)
+scfres = scf.RHF(mol)
+scfres.kernel()
+
+# run ADC calculation using adcc
+state = adcc.run_adc(scfres, method="adc2", n_singlets=2)
+print(state.describe())
+
+# define first SOS term
+Mp_sos_expr = TransitionMoment(f, m_b, n, shifted=True) \
+    * TransitionMoment(n, mup_a, O) / (w_n - w_a)
+Spp_sos_expr = TransitionMoment(f, mup_b, n, shifted=True) \
+    * TransitionMoment(n, mup_a, O) / (w_n - w_a)
+Qpp_sos_expr = TransitionMoment(f, mup_c, n, shifted=True) \
+    * TransitionMoment(n, qp_ab, O) / (w_n - w_a)
+
+# define operator-frequency pairs to be permuted
+Mp_perm_pairs = [(m_b, w_b), (mup_a, w_a)]
+Spp_perm_pairs = [(mup_b, w_b), (mup_a, w_a)]
+Qpp_perm_pairs = [(mup_c, w_b), (qp_ab, w_a)]
+
+tensors = {
+    "Mp": (Mp_sos_expr, Mp_perm_pairs),
+    "Spp": (Spp_sos_expr, Spp_perm_pairs),
+    "Qpp": (Qpp_sos_expr, Qpp_perm_pairs),
+}
+
+# define gauge origin for gauge origin dependent operator integrals
+gauge_origin = "mass_center"
+tpcd_data = {}
+
+# compute two-photon tensors 
+for n_f in range(2):
+    data_state = {}
+    for tensor_name, (sos_expr, perm_pairs) in tensors.items():
+        tensor = compute_tp_tensor(state, sos_expr, n_f, perm_pairs, gauge_origin)
+        if len(tensor.shape) == 3:
+            data_state[tensor_name] = np.einsum("bcd,acd->ab", epsilon, tensor)
+        else:
+            data_state[tensor_name] = tensor
+
+    tpcd_data[f"state_{n_f}"] = data_state
+
+tpcd_data["excitation_energy_uncorrected"] = state.excitation_energy_uncorrected
+
+for ex_state in range(2):
+    R_TP_p = compute_rotatory_strength(tpcd_data, ex_state)
+    print("The two-photon rotatory strength in velocity gauge for excited state "
+          f"{ex_state} is {R_TP_p:.2f} (a.u.).")
+    
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index dd75f4d..c5e633a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -10,7 +10,7 @@ description = "Fun with response function in the ADC framework."
 version = "0.3.0"
 readme = "README.md"
 authors = [
-    { name = "Antonia Papapostolou, Maximilian Scheurer" }
+    { name = "Antonia Papapostolou, Maximilian Scheurer, Friederike Schneider" }
 ]
 license = {file = "LICENSE"}
 # See https://pypi.org/classifiers/
@@ -22,7 +22,7 @@ requires-python = ">=3.7"
 # Declare any run-time dependencies that should be installed with the package.
 dependencies = [
     "adcc>=0.16.0",
-    "respondo>=0.0.5",
+    "respondo>=0.0.6",
     "numpy>=1.14",
     "sympy<=1.13",
     "tqdm",
diff --git a/responsefun/test_property.py b/responsefun/test_property.py
index 2ad7b2d..a4d2cf4 100644
--- a/responsefun/test_property.py
+++ b/responsefun/test_property.py
@@ -3,6 +3,7 @@
 import pytest
 from adcc.Excitation import Excitation
 from adcc.misc import assert_allclose_signfix
+from respondo.mcd import mcd_bterm
 from respondo.polarizability import (
     complex_polarizability,
     real_polarizability,
@@ -16,14 +17,16 @@
     evaluate_property_sos,
     evaluate_property_sos_fast,
 )
-from responsefun.misc import ev2au
+from responsefun.misc import epsilon, ev2au
 from responsefun.SumOverStates import TransitionMoment
 from responsefun.symbols_and_labels import (
     O,
     f,
     gamma,
+    j,
     k,
     m,
+    m_c,
     mu_a,
     mu_b,
     mu_c,
@@ -35,6 +38,7 @@
     w_2,
     w_3,
     w_f,
+    w_j,
     w_k,
     w_m,
     w_n,
@@ -91,6 +95,20 @@ def run_scf(molecule, basis, backend="pyscf"):
         ),
         None,
     ),
+    "mcd_term1": (
+        (
+        TransitionMoment(O, m_c, k) * TransitionMoment(k, mu_b, j, shifted=True)
+        * TransitionMoment(j, mu_a, O) / w_k
+        ),
+        None
+    ),
+    "mcd_term2": (
+        (
+        TransitionMoment(O, mu_b, k) * TransitionMoment(k, m_c, j) * TransitionMoment(j, mu_a, O)
+        / (w_k - w_j)
+        ),
+        None
+    ),
     "beta": (
         (
             TransitionMoment(O, mu_a, n)
@@ -141,7 +159,6 @@ def run_scf(molecule, basis, backend="pyscf"):
     ),
 }
 
-# TODO: add mcd test as soon as gator-program/respondo#15 is merged
 @pytest.mark.parametrize("case", cache.cases)
 class TestIsrAgainstRespondo:
     def test_static_polarizability(self, case):
@@ -253,6 +270,34 @@ def test_tpa_resonant(self, case):
                 tpa, tpa_ref[1], atol=1e-7, err_msg="final_state = {}".format(final_state)
             )
 
+    def test_mcd(self, case):
+        molecule, basis, method = case.split("_")
+        scfres = run_scf(molecule, basis)
+        refstate = adcc.ReferenceState(scfres)
+        state = adcc.run_adc(refstate, method=method, n_singlets=5)
+        mcd_sos_expr1 = SOS_expressions["mcd_term1"][0]
+        mcd_sos_expr2 = SOS_expressions["mcd_term2"][0]
+
+        for ee in state.excitations:
+            gauge_origin = "origin"
+            final_state = ee.index
+            excited_state = Excitation(state, final_state, method)
+            bterm_ref = mcd_bterm(excited_state, gauge_origin=gauge_origin)
+            mcd_tens1 = evaluate_property_isr(
+                state, mcd_sos_expr1, [k],
+                excluded_states=O, excited_state=final_state, gauge_origin=gauge_origin
+            )
+            mcd_tens2 = evaluate_property_isr(
+                state, mcd_sos_expr2, [k],
+                excluded_states=[O,j], excited_state=final_state, gauge_origin=gauge_origin
+            )
+
+            bterm = np.einsum("abc,abc->", epsilon, mcd_tens1 + mcd_tens2)
+
+            np.testing.assert_allclose(
+                bterm, bterm_ref, atol=1e-5, err_msg="final_state = {}".format(final_state)
+            )
+
 
 @pytest.mark.parametrize("case", [case for case in cache.cases if case in cache.data_fulldiag])
 class TestIsrAgainstSos: