diff --git a/qiskit/synthesis/__init__.py b/qiskit/synthesis/__init__.py
index 4d118ddb1fa3..ff7a3b025a8f 100644
--- a/qiskit/synthesis/__init__.py
+++ b/qiskit/synthesis/__init__.py
@@ -46,6 +46,7 @@
    :toctree: ../stubs/
 
     synth_cz_depth_line_mr
+    synth_cx_cz_depth_line_my
 
 Permutation Synthesis
 =====================
@@ -118,7 +119,7 @@
     synth_cnot_count_full_pmh,
     synth_cnot_depth_line_kms,
 )
-from .linear_phase import synth_cz_depth_line_mr, synth_cnot_phase_aam
+from .linear_phase import synth_cz_depth_line_mr, synth_cx_cz_depth_line_my, synth_cnot_phase_aam
 from .clifford import (
     synth_clifford_full,
     synth_clifford_ag,
diff --git a/qiskit/synthesis/clifford/clifford_decompose_layers.py b/qiskit/synthesis/clifford/clifford_decompose_layers.py
index f4f8a1f6c7b3..e859219ba0f2 100644
--- a/qiskit/synthesis/clifford/clifford_decompose_layers.py
+++ b/qiskit/synthesis/clifford/clifford_decompose_layers.py
@@ -23,7 +23,7 @@
     synth_cnot_depth_line_kms,
 )
 from qiskit.synthesis.linear_phase import synth_cz_depth_line_mr
-
+from qiskit.synthesis.linear_phase.cx_cz_depth_lnn import synth_cx_cz_depth_line_my
 from qiskit.synthesis.linear.linear_matrix_utils import (
     calc_inverse_matrix,
     _compute_rank,
@@ -135,10 +135,17 @@ def synth_clifford_layers(
     )
 
     layeredCircuit.append(S2_circ, qubit_list)
-    layeredCircuit.append(CZ2_circ, qubit_list)
 
-    CXinv = CX_circ.copy().inverse()
-    layeredCircuit.append(CXinv, qubit_list)
+    if cx_cz_synth_func is None:
+        layeredCircuit.append(CZ2_circ, qubit_list)
+
+        CXinv = CX_circ.copy().inverse()
+        layeredCircuit.append(CXinv, qubit_list)
+
+    else:
+        # note that CZ2_circ is None and built into the CX_circ when
+        # cx_cz_synth_func is not None
+        layeredCircuit.append(CX_circ, qubit_list)
 
     layeredCircuit.append(H2_circ, qubit_list)
     layeredCircuit.append(S1_circ, qubit_list)
@@ -347,10 +354,17 @@ def _decompose_hadamard_free(
             S2_circ.s(i)
 
     if cx_cz_synth_func is not None:
-        CZ2_circ, CX_circ = cx_cz_synth_func(
-            destabz_update, cliff.destab_x.transpose(), num_qubits=num_qubits
-        )
-        return S2_circ, CZ2_circ, CX_circ
+        # The cx_cz_synth_func takes as input Mx/Mz representing a CX/CZ circuit
+        # and returns the circuit -CZ-CX- implementing them both
+        for i in range(num_qubits):
+            destabz_update[i][i] = 0
+
+        mat_z = destabz_update
+        mat_x = calc_inverse_matrix(destabx.transpose())
+
+        CXCZ_circ = cx_cz_synth_func(mat_x, mat_z)
+
+        return S2_circ, QuantumCircuit(num_qubits), CXCZ_circ
 
     CZ2_circ = cz_synth_func(destabz_update)
 
@@ -399,7 +413,7 @@ def _calc_pauli_diff(cliff, cliff_target):
 def synth_clifford_depth_lnn(cliff):
     """Synthesis of a Clifford into layers for linear-nearest neighbour connectivity.
 
-    The depth of the synthesized n-qubit circuit is bounded by 9*n+4, which is not optimal.
+    The depth of the synthesized n-qubit circuit is bounded by 7*n+2, which is not optimal.
     It should be replaced by a better algorithm that provides depth bounded by 7*n-4 [3].
 
     Args:
@@ -423,6 +437,7 @@ def synth_clifford_depth_lnn(cliff):
         cliff,
         cx_synth_func=synth_cnot_depth_line_kms,
         cz_synth_func=synth_cz_depth_line_mr,
+        cx_cz_synth_func=synth_cx_cz_depth_line_my,
         cz_func_reverse_qubits=True,
     )
     return circ
diff --git a/qiskit/synthesis/linear_phase/__init__.py b/qiskit/synthesis/linear_phase/__init__.py
index 0be6c9c91175..a73f956e5015 100644
--- a/qiskit/synthesis/linear_phase/__init__.py
+++ b/qiskit/synthesis/linear_phase/__init__.py
@@ -13,4 +13,5 @@
 """Module containing cnot-phase circuits"""
 
 from .cz_depth_lnn import synth_cz_depth_line_mr
+from .cx_cz_depth_lnn import synth_cx_cz_depth_line_my
 from .cnot_phase_synth import synth_cnot_phase_aam
diff --git a/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py b/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py
new file mode 100644
index 000000000000..31795bf13aca
--- /dev/null
+++ b/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py
@@ -0,0 +1,262 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2023
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+Given -CZ-CX- transformation (a layer consisting only CNOT gates 
+    followed by a layer consisting only CZ gates)
+Return a depth-5n circuit implementation of the -CZ-CX- transformation over LNN.
+
+Args:
+    mat_z: n*n symmetric binary matrix representing a -CZ- circuit
+    mat_x: n*n invertable binary matrix representing a -CX- transformation
+
+Output:
+    QuantumCircuit: QuantumCircuit object containing a depth-5n circuit to implement -CZ-CX-
+
+References:
+    [1] S. A. Kutin, D. P. Moulton, and L. M. Smithline, "Computation at a distance," 2007.
+    [2] D. Maslov and W. Yang, "CNOT circuits need little help to implement arbitrary 
+        Hadamard-free Clifford transformations they generate," 2022.
+"""
+
+from copy import deepcopy
+import numpy as np
+
+from qiskit.circuit import QuantumCircuit
+from qiskit.synthesis.linear.linear_matrix_utils import calc_inverse_matrix
+from qiskit.synthesis.linear.linear_depth_lnn import _optimize_cx_circ_depth_5n_line
+
+
+def _initialize_phase_schedule(mat_z):
+    """
+    Given a CZ layer (represented as an n*n CZ matrix Mz)
+    Return a scheudle of phase gates implementing Mz in a SWAP-only netwrok
+    (c.f. Alg 1, [2])
+    """
+    n = len(mat_z)
+    phase_schedule = np.zeros((n, n), dtype=int)
+    for i, j in zip(*np.where(mat_z)):
+        if i >= j:
+            continue
+
+        phase_schedule[i, j] = 3
+        phase_schedule[i, i] += 1
+        phase_schedule[j, j] += 1
+
+    return phase_schedule
+
+
+def _shuffle(labels, odd):
+    """
+    Args:
+        labels : a list of indices
+        odd : a boolean indicating whether this layer is odd or even,
+    Shuffle the indices in labels by swapping adjacent elements
+    (c.f. Fig.2, [2])
+    """
+    swapped = [v for p in zip(labels[1::2], labels[::2]) for v in p]
+    return swapped + labels[-1:] if odd else swapped
+
+
+def _make_seq(n):
+    """
+    Given the width of the circuit n,
+    Return the labels of the boxes in order from left to right, top to bottom
+    (c.f. Fig.2, [2])
+    """
+    seq = []
+    wire_labels = list(range(n - 1, -1, -1))
+
+    for i in range(n):
+        wire_labels_new = (
+            _shuffle(wire_labels, n % 2)
+            if i % 2 == 0
+            else wire_labels[0:1] + _shuffle(wire_labels[1:], (n + 1) % 2)
+        )
+        seq += [
+            (min(i), max(i)) for i in zip(wire_labels[::2], wire_labels_new[::2]) if i[0] != i[1]
+        ]
+        wire_labels = wire_labels_new
+
+    return seq
+
+
+def _swap_plus(instructions, seq):
+    """
+    Given CX instructions (c.f. Thm 7.1, [1]) and the labels of all boxes,
+    Return a list of labels of the boxes that is SWAP+ in descending order
+        * Assumes the instruction gives gates in the order from top to bottom,
+          from left to right
+        * SWAP+ is defined in section 3.A. of [2]. Note the northwest
+          diagonalization procedure of [1] consists exactly n layers of boxes,
+          each being either a SWAP or a SWAP+. That is, each northwest
+          diagonalization circuit can be uniquely represented by which of its
+          n(n-1)/2 boxes are SWAP+ and which are SWAP.
+    """
+    instr = deepcopy(instructions)
+    swap_plus = set()
+    for i, j in reversed(seq):
+        cnot_1 = instr.pop()
+        instr.pop()
+
+        if instr == [] or instr[-1] != cnot_1:
+            # Only two CNOTs on same set of controls -> this box is SWAP+
+            swap_plus.add((i, j))
+        else:
+            instr.pop()
+    return swap_plus
+
+
+def _update_phase_schedule(n, phase_schedule, swap_plus):
+    """
+    Given phase_schedule initialized to induce a CZ circuit in SWAP-only network and list of SWAP+ boxes
+    Update phase_schedule for each SWAP+ according to Algorithm 2, [2]
+    """
+    layer_order = list(range(n))[-3::-2] + list(range(n))[-2::-2][::-1]
+    order_comp = np.argsort(layer_order[::-1])
+
+    # Go through each box by descending layer order
+
+    for i in layer_order:
+        for j in range(i + 1, n):
+            if (i, j) not in swap_plus:
+                continue
+            # we need to correct for the effected linear functions:
+
+            # We first correct type 1 and type 2 by switching
+            # the phase applied to c_j and c_i+c_j
+            phase_schedule[j, j], phase_schedule[i, j] = phase_schedule[i, j], phase_schedule[j, j]
+
+            # Then, we go through all the boxes that permutes j BEFORE box(i,j) and update:
+
+            for k in range(n):  # all boxes that permutes j
+                if k in (i, j):
+                    continue
+                if (
+                    order_comp[min(k, j)] < order_comp[i]
+                    and phase_schedule[min(k, j), max(k, j)] % 4 != 0
+                ):
+                    phase = phase_schedule[min(k, j), max(k, j)]
+                    phase_schedule[min(k, j), max(k, j)] = 0
+
+                    # Step 1, apply phase to c_i, c_j, c_k
+                    for l_s in (i, j, k):
+                        phase_schedule[l_s, l_s] = (phase_schedule[l_s, l_s] + phase * 3) % 4
+
+                    # Step 2, apply phase to c_i+ c_j, c_i+c_k, c_j+c_k:
+                    for l1, l2 in [(i, j), (i, k), (j, k)]:
+                        ls = min(l1, l2)
+                        lb = max(l1, l2)
+                        phase_schedule[ls, lb] = (phase_schedule[ls, lb] + phase * 3) % 4
+    return phase_schedule
+
+
+def _apply_phase_to_nw_circuit(n, phase_schedule, seq, swap_plus):
+    """
+    Given
+        Width of the circuit (int n)
+        A CZ circuit, represented by the n*n phase schedule phase_schedule
+        A CX circuit, represented by box-labels (seq) and whether the box is SWAP+ (swap_plus)
+            *   This circuit corresponds to the CX tranformation that tranforms a matrix to
+                a NW matrix (c.f. Prop.7.4, [1])
+            *   SWAP+ is defined in section 3.A. of [2].
+            *   As previously noted, the northwest diagonalization procedure of [1] consists
+                of exactly n layers of boxes, each being either a SWAP or a SWAP+. That is,
+                each northwest diagonalization circuit can be uniquely represented by which
+                of its n(n-1)/2 boxes are SWAP+ and which are SWAP.
+    Return a QuantumCircuit that computes the phase scheudle S inside CX
+    """
+    cir = QuantumCircuit(n)
+
+    wires = list(zip(range(n), range(1, n)))
+    wires = wires[::2] + wires[1::2]
+
+    for i, (j, k) in zip(range(len(seq) - 1, -1, -1), reversed(seq)):
+        w1, w2 = wires[i % (n - 1)]
+
+        p = phase_schedule[j, k]
+
+        if (j, k) not in swap_plus:
+            cir.cnot(w1, w2)
+
+        cir.cnot(w2, w1)
+
+        if p % 4 == 0:
+            pass
+        elif p % 4 == 1:
+            cir.sdg(w2)
+        elif p % 4 == 2:
+            cir.z(w2)
+        else:
+            cir.s(w2)
+
+        cir.cnot(w1, w2)
+
+    for i in range(n):
+        p = phase_schedule[n - 1 - i, n - 1 - i]
+        if p % 4 == 0:
+            continue
+        if p % 4 == 1:
+            cir.sdg(i)
+        elif p % 4 == 2:
+            cir.z(i)
+        else:
+            cir.s(i)
+
+    return cir
+
+
+def synth_cx_cz_depth_line_my(mat_x: np.ndarray, mat_z: np.ndarray):
+    """
+    Joint synthesis of a -CZ-CX- circuit for linear nearest neighbour (LNN) connectivity,
+    with 2-qubit depth at most 5n, based on Maslov and Yang.
+    This method computes the CZ circuit inside the CX circuit via phase gate insertions.
+
+    Args:
+        mat_z : a boolean symmetric matrix representing a CZ circuit.
+            Mz[i][j]=1 represents a CZ(i,j) gate
+
+        mat_x : a boolean invertible matrix representing a CX circuit.
+
+    Return:
+        QuantumCircuit : a circuit implementation of a CX circuit following a CZ circuit,
+        denoted as a -CZ-CX- circuit,in two-qubit depth at most 5n, for LNN connectivity.
+
+    Reference:
+        1. Kutin, S., Moulton, D. P., Smithline, L.,
+           *Computation at a distance*, Chicago J. Theor. Comput. Sci., vol. 2007, (2007),
+           `arXiv:quant-ph/0701194 <https://arxiv.org/abs/quant-ph/0701194>`_
+        2. Dmitri Maslov, Willers Yang, *CNOT circuits need little help to implement arbitrary
+           Hadamard-free Clifford transformations they generate*,
+           `arXiv:2210.16195 <https://arxiv.org/abs/2210.16195>`_.
+    """
+
+    # First, find circuits implementing mat_x by Proposition 7.3 and Proposition 7.4 of [1]
+
+    n = len(mat_x)
+    mat_x = calc_inverse_matrix(mat_x)
+
+    cx_instructions_rows_m2nw, cx_instructions_rows_nw2id = _optimize_cx_circ_depth_5n_line(mat_x)
+
+    # Meanwhile, also build the -CZ- circuit via Phase gate insertions as per Algorithm 2 [2]
+    phase_schedule = _initialize_phase_schedule(mat_z)
+    seq = _make_seq(n)
+    swap_plus = _swap_plus(cx_instructions_rows_nw2id, seq)
+
+    _update_phase_schedule(n, phase_schedule, swap_plus)
+
+    qc = _apply_phase_to_nw_circuit(n, phase_schedule, seq, swap_plus)
+
+    for i, j in reversed(cx_instructions_rows_m2nw):
+        qc.cx(i, j)
+
+    return qc
diff --git a/releasenotes/notes/cx_cz_synthesis-3d5ec98372ce1608.yaml b/releasenotes/notes/cx_cz_synthesis-3d5ec98372ce1608.yaml
new file mode 100644
index 000000000000..2f91f60cc002
--- /dev/null
+++ b/releasenotes/notes/cx_cz_synthesis-3d5ec98372ce1608.yaml
@@ -0,0 +1,32 @@
+---
+features:
+  - |
+    Added a new synthesis algorithm :func:`qiskit.synthesis.linear_phase.synth_cx_cz_depth_line_my`
+    of a CX circuit followed by a CZ circuit for linear nearest neighbor (LNN) connectivity in
+    2-qubit depth of at most 5n using CX and phase gates (S, Sdg or Z). The synthesis algorithm is
+    based on the paper of Maslov and Yang (https://arxiv.org/abs/2210.16195).
+    The algorithm accepts a binary invertible matrix ``mat_x`` representing the CX-circuit,
+    a binary symmetric matrix ``mat_z`` representing the CZ-circuit, and returns a quantum circuit
+    with 2-qubit depth of at most 5n computing the composition of the CX and CZ circuits.
+    The following example illustrates the new functionality::
+
+        import numpy as np
+        from qiskit.synthesis.linear_phase import synth_cx_cz_depth_line_my
+        mat_x = np.array([[0, 1], [1, 1]])
+        mat_z = np.array([[0, 1], [1, 0]])
+        qc = synth_cx_cz_depth_line_my(mat_x, mat_z)
+
+    This algorithm is now used by default in the Clifford synthesis algorithm
+    :func:`qiskit.synthesis.clifford.synth_clifford_depth_lnn` that optimizes 2-qubit depth
+    for LNN connectivity, improving the 2-qubit depth from 9n+4 to 7n+2.
+    The clifford synthesis algorithm can be used as follows::
+
+        from qiskit.quantum_info import random_clifford
+        from qiskit.synthesis import synth_clifford_depth_lnn
+
+        cliff = random_clifford(3)
+        qc = synth_clifford_depth_lnn(cliff)
+
+    The above synthesis can be further improved as described in the paper by Maslov and Yang,
+    using local optimization between 2-qubit layers. This improvement is left for follow-up
+    work.
\ No newline at end of file
diff --git a/test/python/synthesis/test_clifford_decompose_layers.py b/test/python/synthesis/test_clifford_decompose_layers.py
index 79b791141188..cf15f11a1f90 100644
--- a/test/python/synthesis/test_clifford_decompose_layers.py
+++ b/test/python/synthesis/test_clifford_decompose_layers.py
@@ -10,6 +10,7 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
+
 """Tests for Clifford synthesis methods."""
 
 import unittest
@@ -58,11 +59,11 @@ def test_decompose_lnn_depth(self, num_qubits):
         for _ in range(samples):
             cliff = random_clifford(num_qubits, seed=rng)
             circ = synth_clifford_depth_lnn(cliff)
-            # Check that the Clifford circuit 2-qubit depth is bounded by 9*n+4
+            # Check that the Clifford circuit 2-qubit depth is bounded by 7*n+2
             depth2q = (circ.decompose()).depth(
                 filter_function=lambda x: x.operation.num_qubits == 2
             )
-            self.assertTrue(depth2q <= 9 * num_qubits + 4)
+            self.assertTrue(depth2q <= 7 * num_qubits + 2)
             # Check that the Clifford circuit has linear nearest neighbour connectivity
             self.assertTrue(check_lnn_connectivity(circ.decompose()))
             cliff_target = Clifford(circ)
diff --git a/test/python/synthesis/test_cx_cz_synthesis.py b/test/python/synthesis/test_cx_cz_synthesis.py
new file mode 100644
index 000000000000..eec3afe1dfff
--- /dev/null
+++ b/test/python/synthesis/test_cx_cz_synthesis.py
@@ -0,0 +1,84 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Test -CZ-CX- joint synthesis function."""
+
+import unittest
+from test import combine
+import numpy as np
+from ddt import ddt
+from qiskit import QuantumCircuit
+
+from qiskit.quantum_info import Clifford
+
+from qiskit.synthesis.linear_phase.cx_cz_depth_lnn import synth_cx_cz_depth_line_my
+from qiskit.synthesis.linear import (
+    synth_cnot_depth_line_kms,
+    random_invertible_binary_matrix,
+)
+
+from qiskit.synthesis.linear.linear_circuits_utils import check_lnn_connectivity
+
+from qiskit.test import QiskitTestCase
+
+
+@ddt
+class TestCXCZSynth(QiskitTestCase):
+    """Test the linear reversible circuit synthesis functions."""
+
+    @combine(num_qubits=[3, 4, 5, 6, 7, 8, 9, 10])
+    def test_cx_cz_synth_lnn(self, num_qubits):
+        """Test the CXCZ synthesis code for linear nearest neighbour connectivity."""
+        seed = 1234
+        rng = np.random.default_rng(seed)
+        num_gates = 10
+        num_trials = 8
+
+        for _ in range(num_trials):
+            # Generate a random CZ circuit
+            mat_z = np.zeros((num_qubits, num_qubits))
+            cir_z = QuantumCircuit(num_qubits)
+            for _ in range(num_gates):
+                i = rng.integers(num_qubits)
+                j = rng.integers(num_qubits)
+                if i != j:
+                    cir_z.cz(i, j)
+                    if j > i:
+                        mat_z[i][j] = (mat_z[i][j] + 1) % 2
+                    else:
+                        mat_z[j][i] = (mat_z[j][i] + 1) % 2
+
+            # Generate a random CX circuit
+            mat_x = random_invertible_binary_matrix(num_qubits, seed=rng)
+            mat_x = np.array(mat_x, dtype=bool)
+            cir_x = synth_cnot_depth_line_kms(mat_x)
+
+            # Joint Synthesis
+
+            cir_zx_test = QuantumCircuit.compose(cir_z, cir_x)
+
+            cir_zx = synth_cx_cz_depth_line_my(mat_x, mat_z)
+
+            # Check that the output circuit 2-qubit depth is at most 5n
+
+            depth2q = cir_zx.depth(filter_function=lambda x: x.operation.num_qubits == 2)
+            self.assertTrue(depth2q <= 5 * num_qubits)
+
+            # Check that the output circuit has LNN connectivity
+            self.assertTrue(check_lnn_connectivity(cir_zx))
+
+            # Assert that we get the same elements as other methods
+            self.assertEqual(Clifford(cir_zx), Clifford(cir_zx_test))
+
+
+if __name__ == "__main__":
+    unittest.main()