Pluto SDR TX1 & TX2 random phase shift on each measurement

I am currently working with the ADALM-Pluto SDR in its 2R2T configuration and facing an issue regarding phase consistency between two transmitted signals. My test involves transmitting two sinusoidal signals with a predefined phase difference using the PlutoSDR. During it, each transmit channel (TX1 and TX2) is connected to its corresponding receive channel (RX1 and RX2).

Each time I get random phase offset, for example, two consecutive tests (expected phase shift 30 degrees):

With loopback enabled ("sdr.loopback = 1") I get the expected result:

I understand that RX and TX have separate LOs, but this should not affect the phase difference between the TX signals.

What am I doing wrong? Does such a task require some specific calibration?

Code I am using for the test:

SDR setup:

# Configuration parameters
sample_rate = 20e6 # must be <=30.72 MHz if both channels are enabled
rf_bandwidth = 1e6  # RF bandwidth 
carrier_freq = 400e6  # Frequency of the tone
tx_gain0 = -30  # TX gain in dB
tx_gain1 = -30  # TX gain in dB
num_samples = 2**18  # Number of samples in the buffer
rx_mode = "manual"  # can be "manual" or "slow_attack"
rx_gain0 = 0
rx_gain1 = 0

def initSDR():
    # Connect to Pluto SDR
    sdr = adi.ad9361(uri=PLUTO_IP)

    # Configure SDR properties
    sdr.gain_control_mode = rx_mode
    sdr.rx_rf_bandwidth = int(rf_bandwidth)  # RX bandwidth
    sdr.sample_rate = int(sample_rate)  # RX sample rate 
    sdr.rx_lo = int(carrier_freq)  # RX local oscillator frequency
    sdr.rx_hardwaregain_chan0 = int(rx_gain0)   # RX gain for channel 0
    sdr.rx_hardwaregain_chan1 = int(rx_gain1)   # RX gain for channel 1
    sdr.rx_buffer_size = int(num_samples)   # Set RX buffer size
    sdr.rx_enabled_channels = [0, 1]  # Enable both channels
    
    # Configure TX properties
    sdr.tx_rf_bandwidth = int(rf_bandwidth)  # TX bandwidth
    sdr.tx_lo = int(carrier_freq)  # TX local oscillator frequency
    sdr.tx_hardwaregain_chan0 = tx_gain0  # TX gain for channel 0
    sdr.tx_hardwaregain_chan1 = tx_gain1  # TX gain for channel 1
    sdr.tx_cyclic_buffer = True  # Enable cyclic buffer
    sdr.tx_buffer_size = int(num_samples)  # Set TX buffer size
    sdr.tx_enabled_channels = [0, 1]  # Enable both channels

    sdr.loopback = 0    # Enable loopback mode

    return sdr

Signal creation:

def create_signal(phase_shift=0):
    fs = int(sdr.sample_rate)                   # Sampling rate
    # int(fc / (fs / N)) * (fs / N)             # is used to ensure that fc is a multiple of the frequency resolution
    fc0 = int(frequency_of_signal)              # Frequency of the first signal
    fc1 = int(frequency_of_signal)              # Frequency of the second signal
    N = 2**16                                   # Number of samples
    ts = 1 / float(fs)                          # Time step
    t = np.arange(0, N * ts, ts)                # Time vector
    i0 = np.cos(2 * np.pi * t * fc0) * 2 ** 14  # In-phase component of the first signal
    q0 = np.sin(2 * np.pi * t * fc0) * 2 ** 14  # Quadrature component of the first signal
    i1 = np.cos(2 * np.pi * t * fc1 + phase_shift*(np.pi/180)) * 2 ** 14
    q1 = np.sin(2 * np.pi * t * fc1 + phase_shift*(np.pi/180)) * 2 ** 14
    iq0 = i0 + 1j * q0
    iq1 = i1 + 1j * q1
    
    return [iq0, iq1]

Visualization part:

# Collect data from the receiver
data = sdr.rx()
Rx_0 = data[0]  # Data from the first RX channel
Rx_1 = data[1]  # Data from the second RX channel

fs = int(sdr.sample_rate)  # Sampling rate
t = np.arange(len(Rx_0)) / fs  # Time vector

# Define the frequency for time range calculation
frequency_of_waveform = frequency_of_signal  # Frequency of the transmitted signal
period_of_waveform = 1 / frequency_of_waveform
number_of_periods_to_display = 5  # Displaying 5 periods

# Find the time range that covers the specified number of periods
time_range_to_display = period_of_waveform * number_of_periods_to_display
indices_to_display = t < time_range_to_display  # Correcting the range

# Compute phases of both signals
phase_0 = np.angle(Rx_0)
phase_1 = np.angle(Rx_1)

# Calculate phase difference, wrap it and convert to degrees
phase_diff = np.unwrap(phase_1 - phase_0) * (180 / np.pi)
phase_diff = np.mod(phase_diff, 360)  # Ensure phase difference is within 0-360 degrees
mean_phase_diff = np.mean(phase_diff[indices_to_display])

# Plotting both received signals within the correct range
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(t[indices_to_display], np.real(Rx_0)[indices_to_display], label='Channel 0 - TX 0 RX')
plt.plot(t[indices_to_display], np.real(Rx_1)[indices_to_display], label='Channel 1 - TX 1 RX', linestyle='--')
plt.title('Received Signals to Compare Phase Shift')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)

# Plotting phase difference in degrees
plt.subplot(2, 1, 2)
plt.plot(t[indices_to_display], phase_diff[indices_to_display], label='Phase Difference', color='green')
plt.axhline(y=mean_phase_diff, color='r', linestyle='-', label=f'Mean Phase Difference: {mean_phase_diff:.2f}°')
plt.title('Phase Difference Over Time (Degrees)')
plt.xlabel('Time (s)')
plt.ylabel('Phase Difference (Degrees)')
plt.legend()
plt.grid(True)