Merge additional axes in contour

Question

I used the RK4 method to calculate the radiated spectrum for a two- level system in quantum mechanics in matplotlib. I also want to plot the signal along the vertical axis together in the same figure, but my code doesn’t implement it. The code is provided below (I also attached the original figure (Fig. (a)). In general, some clear explanation to plot another x,y-axis as in the first attached figure. Many thanks to all in advance. [enter image description here] [enter image description here] : https://i.stack.imgur.com/lYLCs.jpg. enter image description here The goals of this code are: 1) Numerical calculation of the Bloch equations and 2) Calculated the radiated spectra from a two-level system. Numerical calculation was performed using RK4 method. The Fast Fourier Transform (FFT) was also been used.

    import numpy as np
    from scipy.integrate import odeint
    import matplotlib.pyplot as plot
    import math
    from scipy.fftpack import*
    plot.rcParams.update({'font.size': 16})

    # Parameters
    n= 30                                # Number of optical cycles.
    Om0 = 1.5/27.2114                         # Carrier frequency.
    Om = Om0                                  # Transition frequency.
    dt = 0.1                                  # Time step.
    E0 = 1#np.sqrt(I/I0)                      # Electric peak 
    t = np.arange(0,2*np.pi*n/Om0,0.01)       # Time in x-axis data.
    d = 1.0                                   # Dipole moment.
    h_par = 1.0            # The reduced Planck constant (=1 in a.u.)
    E = E0 * np.cos(Om0*t)                    # Electric pulse.
    OmR = d*E/h_par                        # Rabi Flopping frequency.
    N = np.size(t)                         # Total number of samples.
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ## Numerical solutions using RK4 method.                                      
    # Define Bloch system                      
    def f1(t,u,v,w):
        return Om * v
    def f2(t,u,v,w):
        return -Om * u -2 * (E0*np.cos(Om0*t)) * w
    def f3(t,u,v,w):
        return 2 * (E0*np.cos(Om0*t)) * v
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    def integrate(u0,v0,w0):                                           
        u = np.zeros(N)
        v = np.zeros(N)            # Filled with arrays with zeros.
        w = np.zeros(N)

        u[0] = u0                    # Construct initial conditions.
        v[0] = v0
        w[0] = w0
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        # The fourth-order Runge-Kutta coefficients (RK4)
        for i in range(N-1):        
             k1u = dt * f1(t[i], u[i], v[i], w[i])
             k1v = dt * f2(t[i], u[i], v[i], w[i])
             k1w = dt * f3(t[i], u[i], v[i], w[i])            
             k2u = dt * f1(t[i] +  0.5 * dt, u[i] + 0.5 * k1u, v[i] + 0.5 * k1v, w[i] + 0.5 * k1w)
             k2v = dt * f2(t[i] +  0.5 * dt, u[i] + 0.5 * k1u, v[i] + 0.5 * k1v, w[i] + 0.5 * k1w)
             k2w = dt * f3(t[i] +  0.5 * dt, u[i] + 0.5 * k1u, v[i] + 0.5 * k1v, w[i] + 0.5 * k1w)        
             k3u = dt * f1(t[i] +  0.5 * dt, u[i] + 0.5 * k2u, v[i] + 0.5 * k2v, w[i] + 0.5 * k2w)
             k3v = dt * f2(t[i] +  0.5 * dt, u[i] + 0.5 * k2u, v[i] + 0.5 * k2v, w[i] + 0.5 * k2w)
             k3w = dt * f3(t[i] +  0.5 * dt, u[i] + 0.5 * k2u, v[i] + 0.5 * k2v, w[i] + 0.5 * k2w)            
             k4u = dt * f1(t[i] + dt, u[i] + k3u, v[i] + k3v,w[i] + k3w)
             k4v = dt * f2(t[i] + dt, u[i] + k3u, v[i] + k3v,w[i] + k3w)
             k4w = dt * f3(t[i] + dt, u[i] + k3u, v[i] + k3v,w[i] + k3w)
        
             u[i + 1] = u[i] + (k1u + 2 * k2u + 2 * k3u + k4u)/6
             v[i + 1] = v[i] + (k1v + 2 * k2v + 2 * k3v + k4v)/6
             w[i + 1] = w[i] + (k1w + 2 * k2w + 2 * k3w + k4w)/6            
        return u,v,w 

        # Initial conditions
        u1,v1, w1 = integrate(0.0,0.0,-1.0) 

        ## Fast Fourier Transform (FFT).
        # Frequency domain.
        dom = 2*np.pi/(t[-1]-t[0])                  # Frequency step.
        om = np.arange(-N/2, N/2)*dom                # Generate 
        frequency axis data.
        FFT_u1 = fftshift(fft(u1))*dt/np.sqrt(2*np.pi) 
        S = np.abs(FFT_u1)
        #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        P = np.log10((S)**2)  # Power spectrum of a two-level system.
        #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        # Plotting results
        font = {'family': 'serif',
                'color':  'darkred',
                'weight': 'normal',
                'size': 16,
                }
    
        fig, ax = plot.subplots(1, 1, figsize =(8,6),
                                dpi=300) # Resolution. 
        #ax.plot(om/Om0,P, linestyle='-', c='black')
        #ax.set_xlabel('$\omega/\omega_0$', fontdict=font)
        #ax.set_ylabel(r"$log\mid P \mid^2$", fontdict=font)
        ax.set_xlim(0.0,20.0)
        #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        # Contour plot        
        y = np.linspace(0,20,1000)
        P2 = np.log10((S)**2)*E
        cp = ax.contour(om/Om0,y/Om0, P2)
        fig.colorbar(cp)                   # Add a colorbar to a plot
        fig.suptitle('Numerical Solutions.', fontdict=font)
        plot.tight_layout()
        plot.show()