parosc.xmds

<?xml version="1.0"?>
<!--Parametric oscilator simulation using nu-P (squeezed +P) function-->
<!--Equations have pump adiabatically eliminated, to give 2-photon loss-->
<!--Noise correlations are very complicated, which slows things down-->
<!--Equations are also not very stable (no gauges yet)-->

<simulation>

  <name>parosc</name>

  <author>Joel Corney</author>
  <description>
    Parametric oscilator simulation using nu-P (squeezed +P) function
    Equations have pump adiabatically eliminated, to give 2-photon loss
    Noise correlations are very complicated, which slows things down
    Equations are also not very stable (no gauges yet)
  </description>

  <prop_dim>t</prop_dim>
  <error_check>yes</error_check>
  <stochastic>yes</stochastic>
  <paths>1000</paths>
  <use_mpi>no</use_mpi>
  <seed>1 2</seed>
  <noises>16</noises>

  /***********/
  <globals>
  /***********/
  <![CDATA[
    const double Pi = 2.0*asin(1.0);

    /*Simulation Parameters*/
    const double kappa1 = 1;    /*Damping of alpha-mode*/
    const double kappa2 = 500;  /*Damping of gamma-mode*/
    const double g = 10;        /*Parametric coupling*/
    const double epsilon = 75;   /*Driving*/

    /*Scaled parameters*/
    const double lambda = g*epsilon/kappa2;
    const double h = g/sqrt(2*kappa2);
  ]]>
  </globals>

  /*************/
  <field>
  /*************/
    <name>main</name>
    <samples>1</samples>
 
    /***MAIN***/
    <vector>
      <name>main</name>
      <type>complex</type>
      <components>alpha1 alpha2 nu1 nu2 mu</components>
      <fourier_space>no</fourier_space>    
      <![CDATA[
        /* NB: We take the '2' variables to be the conjugates */
        alpha1 = rcomplex(sqrt(0.0),0.0); //rcomplex(sqrt((lambda- kappa1)/(h*h)),0.0);
        alpha2 = rcomplex(sqrt(0.0),0.0); //rcomplex(sqrt((lambda- kappa1)/(h*h)),0.0);    
        nu1 = rcomplex(0.0,0.0);
        nu2 = rcomplex(0.0,0.0); 
        mu = c_sqrt(1+nu1*nu2);
      ]]>
    </vector>        
  </field>

  /******************/
  <sequence>
  /******************/
    <integrate>
      <algorithm>SIIP</algorithm>
      <interval>2</interval>
      <lattice>5000</lattice>
      <samples>200</samples>
      <iterations>4</iterations>
      <vectors>main</vectors>
      <![CDATA[
      
        complex nunu = nu1*nu2;
        
        complex w2 = rcomplex(n_2,n_11)/sqrt(2.0);      /*Complex noise*/
        complex w3 = rcomplex(n_3,n_12)/sqrt(2.0);      /*Complex noise*/
        complex w4 = rcomplex(n_4,n_13)/sqrt(2.0);      /*Complex noise*/
        complex w6 = rcomplex(n_6,n_14)/sqrt(2.0);      /*Complex noise*/
        complex w7 = rcomplex(n_7,n_15)/sqrt(2.0);      /*Complex noise*/
        complex w9 = rcomplex(n_9,n_16)/sqrt(2.0);      /*Complex noise*/

        complex s1 = c_sqrt(2*mu*nu1*(kappa1 + 2*h*h*alpha1*alpha2));           /*Noise Correlations*/
        complex s5 =  c_sqrt(2*mu*nu2*(kappa1 + 2*h*h*alpha1*alpha2));          /*Noise Correlations*/
        complex s2 =  c_sqrt(-2*nu1*nu2*(kappa1 + 2*h*h*alpha1*alpha2));                /*Noise Correlations*/

        complex s8 =  c_sqrt(-2*h*h*nu1*nu1*(nu1*nu1*nu2*nu2 + 4*nu1*nu2 +2)/(mu*mu));          /*Noise Correlations*/
        complex s10 =  c_sqrt(-2*h*h*nu2*nu2*(nu1*nu1*nu2*nu2 + 4*nu1*nu2 +2)/(mu*mu));         /*Noise Correlations*/
        complex s9 =  c_sqrt(-2*h*h*nu1*nu1*nu1*nu2*nu2*nu2/(mu*mu));           /*Noise Correlations*/

        complex s3 =  c_sqrt(-h*h*nu1*(alpha1*(3*nu1*nu1*nu2*nu2+5*nu1*nu2 +2) - 2*alpha2*mu*nu1)/(mu*mu));             /*Noise Correlations*/
        complex s7 =  c_sqrt(-h*h*nu2*(alpha2*(3*nu1*nu1*nu2*nu2+5*nu1*nu2 +2) - 2*alpha1*mu*nu2)/(mu*mu));             /*Noise Correlations*/
        complex s4 =  c_sqrt(h*h*nu1*nu2*(2*alpha2*(2*nu1*nu2 +1) + alpha1*mu*nu2)/mu);         /*Noise Correlations*/
        complex s6 =  c_sqrt(h*h*nu1*nu2*(2*alpha1*(2*nu1*nu2 +1) + alpha2*mu*nu1)/mu);         /*Noise Correlations*/

        complex nu1_deriv = h*h*(((3*nunu+2)*alpha1*alpha1 - nu1*nu1*alpha2*alpha2)/(2*mu) + nu1*(2*mu*mu*mu*mu + 1)/(2*mu*mu))
                            + lambda*(nu1*nu1 - 3*nunu - 2)/(2*mu)
                             + s8*n_8 + s9*w9 + s3*~w3 + s6*~w6;
        complex nu2_deriv = h*h*(((3*nunu+2)*alpha2*alpha2 - nu2*nu2*alpha1*alpha1)/(2*mu) + nu2*(2*mu*mu*mu*mu + 1)/(2*mu*mu))
                            + lambda*(nu2*nu2 - 3*nunu - 2)/(2*mu)
                             + s10*n_10 + s9*~w9 + s4*~w4 + s7*~w7;

        dalpha1_dt = (h*h/2 - kappa1)*alpha1 + (lambda - h*h*(alpha1*alpha1+3*mu*nu1/2))*alpha2
                      + s1*n_1 + s2*w2 + s3*w3 + s4*w4;
        dalpha2_dt = (h*h/2 - kappa1)*alpha2 + (lambda - h*h*(alpha2*alpha2+3*mu*nu2/2))*alpha1
                      + s5*n_5+ s2*~w2 + s6*w6 + s7*w7;
                   
        dnu1_dt =  nu1_deriv;
        dnu2_dt = nu2_deriv;
        dmu_dt = (nu2*nu1_deriv + nu1*nu2_deriv)/(2*mu);
      ]]>
    </integrate>
  </sequence>

  /**************/
  <output>
    <filename>parosc.xsil</filename>
    <group>
      <sampling>
       <vectors>main</vectors>
        <moments>mumu nunu nunui N_alpha N_alphai alpha_x alpha_y alpha_x2</moments>
        <![CDATA[
          mumu = mu*mu;
          nunu = nu1*nu2;
          N_alpha = alpha1*alpha2 + nunu;
          nunui = imag(nu1*nu2);
          N_alphai = imag(alpha1*alpha2 +nu1*nu2);
          alpha_x = (alpha1 + alpha2)/2;
          alpha_y = (alpha1 - alpha2)/(2*i);
          alpha_x2 = (alpha1 + alpha2)*(alpha1 + alpha2)/4.0 + (1.0 + 2*nunu - mu*(nu1+nu2))/4.0;
        ]]>
      </sampling>
    </group>
   </output>
</simulation>