User:Deciwill

1. pragma rtGlobals=1 // Use modern global access method.

// Optomechanical Theory // V 1.1 D Brooks // optomechanical theory rewritten to work with global variables, fit-compatible functions, and a front panel. // Since math and pi cannot be entered into panel, we will work in cyclic units. // 0 phase shift is antinode. node is at pi*0.5. // v 1.3 N Brahms // Modified for 9-well averaging // v5 N Brahms // Integrated with version 5 code, arbitrary well averaging, some cleanup.

function omFunctionList() // PANELS and INITIALIZATION // makeOmPanel // omInit // CALCULATE OM PARAMETERS // omRecalc // calcAParameters // myprotofunc // avgWell // getContrast // getDN // getKOpt // BroadenSlope // getStaticK // getStaticOsc // getOscFreq // getGQed // getZHo // getLambDicke // getGranularity // getHeating // gOpt // FITTING FUNCTIONS // FreqVsDelta_fit // FreqVsPhi_fit // fano // logFano // BreitWigner // BreitWigner2 // omFitMaster // SSresponse end

Menu "Macros" "Open Optomechanics Panel", makeOmPanel() end

Window makeOmPanel(): Panel DoWindow /k omPanel NewPanel /w=(370,320,370+520,320+230) /N=omPanel /k=1 as "Optomechanics Panel" //Optomechanic variables variable/g OptoMechInit //Initializes OptoMechanic variables only first time panel is called if(!OptoMechInit) omInit("") endif

//Other parameters //variable/g sLa0, sLa1p, sLa1n, sLa2p,sLa2n,sLn // linear sensitivity for five well weightings

make/o/t/n=(numpnts(omA)) omAList omAList = "A["+num2str(p)+"] = "+num2str(omA[p]) // variable i // for (i=0; i<numpnts(omA); i+=1) // omAList[p] = num2str(omA[p]); // endfor

//column 1 variable c=10; variable r=5; button resetparams, pos={c,r}, size={150,30}, proc=omInit, title="Click to Reset\r Parameters to Default"; r+=40 setvariable omG0, pos={c,r}, size={150,15}, proc=omRecalc, title="g_0 [Hz] ", value=omG0, bodywidth=70, format="%4.3g", disable=2; r+=20 setvariable omK, pos={c,r}, size={150,15}, proc=omRecalc, title="kappa [Hz] ", value=omK, bodywidth=70, format="%4.3g", disable=2; r+=30 setvariable omDCa, pos={c,r}, size={150,15}, proc=omRecalc, title="Delta_ca [Hz] ", value=omDCa, bodywidth=70, format="%4.3g"; r+=20 setvariable omDelta, pos={c,r}, size={150,15}, proc=omRecalc, title="Delta [Hz] ", value=omDelta, bodywidth=70, format="%4.3g"; r+=20 setvariable omNBar, pos={c,r}, size={150,15}, proc=omRecalc, title="Num photons ", value=omNBar, bodywidth=70; r+=20 setvariable omNPhonon, pos={c,r}, size={150,15}, proc=omRecalc, title="Num phonons ", value=omNPhonon, bodywidth=70; r+=20 setvariable omNA, pos={c,r}, size={150,15}, proc=omRecalc, title="Num atoms ", value=omNA, bodywidth=70; r+=20 setvariable omf0,pos={c,r}, size={150,15}, proc=omRecalc, title="f 0 [Hz] ", value=omf0, bodywidth=70, format="%4.3g"; r+=20

//column 2 c+=170; r=5; setvariable omPhi, pos={c,r}, size={150,15}, proc=omRecalc, title="Phi [pi] ", value=omPhi, bodywidth=70; r+=20 setvariable omNoise, pos={c,r}, size={150,15}, title="Cavity noise [Hz] ", value=omDeltaNoise, bodywidth=70, format="%4.3g"; r+=20; setvariable omASize, pos={c,r}, size={150,15}, proc=omRecalc, title="Dist. width [a]", value=omASize, bodywidth=70; r+=20 setvariable omNumA, pos={c,r}, size={150,15}, proc=omRecalc, title="# avg. coeffs. ", value=omNumA, bodywidth=70; r+=30 valdisplay omContrast,pos={c,r}, size={140,15}, title="Contrast: ", value=omContrast, bodywidth=50, format="%4.3f"; r+=20 valdisplay omDN, pos={c,r}, size={140,15}, title="Delta_N [MHz]: ", value=omDN, bodywidth=50, format="%4.1f"; r+=25 listbox omAList, pos={c+20,r}, size={130,72}, listWave = omAList; // valdisplay omA0,pos={c+5,r},size={85,15},title="a0 ",value=omA0;r+=20 // valdisplay omA1,pos={c+5,r},size={85,15},title="a1 ",value=omA1;r+=20 // valdisplay omA2,pos={c+5,r},size={85,15},title="a2 ",value=omA2;r+=20

//column 3 c+=180; r=5 titlebox omFreqTitle, pos={c,r}, size={150,15}, frame=0, title="Oscillation frequencies (Hz):"; r+=20 valdisplay omOmegaOsc,pos={c,r}, size={150,15}, title="Total ", value=omOmegaOsc, bodywidth=60, format="%4.3g"; r+=20 valdisplay omOmegaStatic,pos={c,r}, size={150,15}, title="Static ", value=omOmegaStatic,bodywidth=60, format="%4.3g";r+=20 valdisplay omOmegaShift,pos={c,r}, size={150,15}, title="Shift ", value=omOmegaShift, bodywidth=60, format="%4.3g"; r+=30

valdisplay omGQed, pos={c,r}, size={150,15}, title="cQED coupling (Hz) ", value=omGQED, bodywidth=60, format="%4.3g"; r+=20 valdisplay omGOm, pos={c,r}, size={150,15}, title="OM coupling (Hz) ", value=omGOm, bodywidth=60, format="%4.3g"; r+=20 valdisplay omZSize, pos={c,r}, size={150,15}, title="Z size (nm) ", value=omZSize, bodywidth=60, format="%4.1f"; r+=20 valdisplay omEpsilon, pos={c,r}, size={150,15}, title="Granularity param ", value=omEpsilon, bodywidth=60, format="%4.3f"; r+=20 valdisplay omHeating, pos={c,r}, size={150,15}, title="Heating (quanta/s) ", value=omHeating, bodywidth=60, format="%4.3g"; r+=20

omRecalc("",0,"","") end

//return to omFunctionList //-----------------------------------------------------------------------------------------------

//----------------------------------------------------------------------------------------------- function omInit(ctrlname) : ButtonControl String ctrlName variable/g omG0 = 13.1e6 // per-atom on resonance coupling strength variable/g omK = 1.82e6 // cavity linewidth via ringdown measurement variable/g omPhi = 1/4 // probe-trap phase shift at atom well, 0 is antinode variable/g omNBar = 3.5 // mean photon number, averaged over the sidelock bandwidth variable/g omNA = 4000 // number of atoms variable/g omNPhonon = 0.15 // number of thermal COM phonons variable/g omDCa = -40e9 // detuning of bare cavity resonance from atomic resonance (red is negative) variable/g omMass = 1.44466881e-25 // per-atom mass of Rb, m=87 amu variable/g omKP = 2*pi/(780e-9) // probe wavenumber variable/g omKT = 2*pi/(850e-9) // trap wavenumber variable/g omf0 = 57.7e3 // trap frequency variable/g omfs = 140e3 // shifted frequency variable/g omDelta = -2.05e6 // detuning of probe from averaged with-atom cavity resonance (red is negative) variable/g omASize = .88 // Size of cloud width in units of trap wells (425 nm) variable/g omNumA = 3 // Number of wells to include in well-weighting (must be integer > 1) make/o/n=(omNumA) omA // Site averaging coefficients make/o/n=(2*omNumA-1) omL // Linear sensitivity coefficients //variable/g omA0,omA1,omA2,omA3,omA4 variable/g omDPhi = 0.27 /pi // phi spacing for multi-well averaging; There are 11.64 wells per phase of pi variable/g omDeltaNoise = .6e6 // Standard Deviation of the gaussian technical noise that broadens cavity linewidths and slopes.

variable/g omContrast variable/g omOmegaOsc variable/g omOmegaStatic variable/g omOmegaShift

variable/g omGQed // Effective cavity QED coupling parameter, = sqrt(N)*g0^2/dca variable/g omGOm // Effective optomechanical coupling parameter, = sqrt(N)*g0^2/dca*kp*zho variable/g omZSize // Size within well (nm) variable/g omEpsilon // Granularity parameter variable/g omHeating // Per-atom heating rate (quanta/s) variable/g omDN // Delta N (MHz)

omRecalc("",0,"","") nvar OptoMechInit OptoMechInit=1 end

//return to omFunctionList //------------------------------------------------------------------------------------------------ // function omRecalc(ctrlName, varNum,varStr, varName) // recalculates omContrast and linear sensitivity parameters when a variable is changed on the front panel // can also be called from functions, useful during fit routines //------------------------------------------------------------------------------------------------ function omRecalc(ctrlName, varNum,varStr, varName) :SetVariableControl String ctrlName Variable varNum // value of variable as number String varStr // value of variable as string String varName // name of variable

calcAParameters() nvar omContrast, omASize, omPhi omContrast=getContrast()

nvar omOmegaOsc, omOmegaStatic, omOmegaShift, omDelta,omNA,omf0 nvar omGQed, omGOm, omZSize, omEpsilon, omHeating, omKP, omNPhonon, omDN

omOmegaOsc = getOscfreq(omDelta,omNA,omf0,omASize,omPhi,initialized=1) omOmegaStatic = getStaticOsc(omPhi) omOmegaShift=Sign(omOmegaOsc^2-omOmegaStatic^2)*(Abs(omOmegaOsc^2-omOmegaStatic^2))^.5

omGQed = getGQed(); omGOm = omKP*getZHO()*omGQed*sqrt(omNA); // omLD = getLambDicke(); omZSize = (2*omNPhonon+1)*getZHO()*1e9; omDN = avgWell(getDN)/1e6; omEpsilon = getGranularity(); omHeating = avgWell(getHeating);

DoUpdate /W=omPanel end

//return to omFunctionList //------------------------------------------------------------------------------------------------ // function calcAParameters //------------------------------------------------------------------------------------------------ function calcAParameters() nvar omASize,omNumA,omPhi,omDPhi

// Need to average over at least one well if (mod(omNumA,1)!=0 || omNumA < 1) DoAlert 0, "Number of wells must be an integer >= 1" wave omA omNumA = numpnts(omA) return -1 endif

// A size must be positive if (omASize < 0) DoAlert 0, "A size must be positive" wave omA; omA = NaN; make/o/t/n=(numpnts(omA)) omAList; omAList = "A["+num2str(p)+"] = NaN"; return -2 endif

make/o/n=(omNumA) omA

// Calculate the well-weighting parameters if (omASize==0) omA = 0; omA[0] = 1; elseif (omNumA>1) variable iA, aTotal = 0 for(iA=0; iA<omNumA-1; iA+=1) omA[iA] = .5*Erf((iA+.5)/omASize/sqrt(2)) - aTotal; aTotal += omA[iA]; endfor omA[iA] = .5*Erfc((iA-.5)/omASize/sqrt(2)); // Last wells omA[0] *= 2; // Correct for only having one center well else omA[0] = 1; endif

// Make labels for the list box make/o/t/n=(numpnts(omA)) omAList omAList = "A["+num2str(p)+"] = "+num2str(omA[p])

// Calculate linear sensitivity make/o/n=(2*omNumA-1) omL omL[omNumA-1] = omA[0]*sin(2*pi*omPhi) // Handle center well for(iA=1; iA<numpnts(omA); iA+=1) omL[omNumA-iA-1] = omA[iA]*sin(2*pi*(omPhi-iA*omDPhi)) omL[omNumA+iA-1] = omA[iA]*sin(2*pi*(omPhi+iA*omDPhi)) endfor

end

//return to omFunctionList //------------------------------------------------------------------------------------------------ // function myprotofunc(localPhi) // serves as a template for functions which are to be passed to another function using FUNCREF // takes a phase 0 to 0.5 as an input // returns some calculation based on that omPhi //------------------------------------------------------------------------------------------------ function myprotofunc(localPhi) variable localPhi nvar omDPhi return localPhi+omDPhi end

//------------------------------------------------------------------------------------------------ // function avgWell(fin,[localPhi]) // Averages input function, fin, over omNumA wells weighted by atom density in those wells // if a omPhi isn't specified, uses the global omPhi value //------------------------------------------------------------------------------------------------ function avgWell(fin,[phi]) FUNCREF myprotofunc fin variable phi nvar omDPhi,omPhi wave omA

//Default Parameters if(ParamisDefault(phi)) phi=omPhi endif

variable avg = omA[0]*fin(phi) // Center well variable iA for (iA=1; iA<numpnts(omA); iA+=1) // All other wells avg += omA[iA]*(fin(phi+iA*omDPhi)+fin(phi-iA*omDPhi)) endfor return avg end //return to omFunctionList

//------------------------------------------------------------------------------------------------ // getContrast returns the 5 well averaged omContrast //------------------------------------------------------------------------------------------------ function getContrast()

variable DNmax, DNmin

DNmax = avgWell(getDN,phi=0) DNmin = avgWell(getDN,phi=0.5)

return (DNmax-DNmin)/(DNmax+DNmin) end //------------------------------------------------------------------------------------------------ // getDN calculates the atom-cavity shifted resonance for the atoms at phase localphi //------------------------------------------------------------------------------------------------ function getDN(phi) variable phi variable DNphipos nvar omNA,omG0,omDCa,hbar,omKP,omMass,omf0,omNPhonon

DNphipos=(omNA*omG0^2/omDCa)*.5*( 1 + exp(-getLambDicke()^2) * cos(2*pi*phi) )

return DNphipos end //return to omFunctionList

//------------------------------------------------------------------------------------------------ // getKOpt calculates the oscillation frequency of the atoms located at local omPhi //------------------------------------------------------------------------------------------------ function getKOpt(localPhi) variable localPhi variable statick, dynamickprefactor, dynamick nvar omDelta,omNA,omf0,omASize,omPhi nvar hbar, omG0, omKP, omDCa, omNBar, omDPhi,omK,omMass,omDeltaNoise

statick = -2 * hbar * (2*pi*omG0)^2 * omKP^2 * omNBar / (2*pi*omDCa)*cos(2*pi*LocalPhi) dynamickprefactor = 2 * hbar * omKP^2 * omNA * (2*pi*getGQed())^2 * omNBar* sin(2*pi*LocalPhi)^2 dynamick=dynamickprefactor*Integrate1d(BroadenSlope,(-5*omDeltaNoise),(5*omDeltaNoise)) return (statick+dynamick) end

//------------------------------------------------------------------------------------------------ //BroadenSlope is used to calculate noise-broaden slope to be used for dynamic shifts //------------------------------------------------------------------------------------------------ function BroadenSlope(inx) variable inx nvar omDelta, omK, omDeltaNoise

variable slopefunc=((2*pi*(omDelta+inx)) / ((2*pi*omK)^2 + (2*pi*(omDelta+inx))^2) return 1/sqrt(2*pi)/(omDeltaNoise)*Exp(-(inx^2)/(2*omDeltaNoise^2))*slopefunc end //return to omFunctionList

//------------------------------------------------------------------------------------------------ // getStaticK calculates the static shifts only //------------------------------------------------------------------------------------------------ function getStaticK(localPhi) variable localPhi variable statick, dynamick nvar omDelta,omNA,omf0,omASize,omPhi nvar hbar, omG0, omKP, omDCa, omNBar, omDPhi,omK,omMass

statick = -2 * hbar * (2*pi*omG0)^2 * omKP^2 * omNBar / (2*pi*omDCa)*cos(2*pi*LocalPhi) return statick end

//------------------------------------------------------------------------------------------------ // getStaticOsc calculates the oscillation frequency from static shifts only //------------------------------------------------------------------------------------------------ function getStaticOsc(localPhi) variable localPhi nvar omf0,omMass

variable f=sqrt((2*pi*omf0)^2+(avgWell(getStaticK,phi=localPhi))/omMass)/2/pi return f end //return to omFunctionList

//----------------------------------------------------------------------------------------------------- //function getOscFreq(LocDelta,LocNa,LocWt,LocAsize,LocPhi,[initialized]) // Returns the cyclic frequency of the oscillator, // as a function of omDelta, # of atoms, trap frequency, distribution of atoms, and omPhi. // Returns 5-well sensitivity averaged results. // // To be used by fit function, it needs to take local variables, which are then written over the global variables // initialized = 0 by default. use initialized =1 when called from omRecalc //----------------------------------------------------------------------------------------------------- function getOscFreq(LocDelta,LocNa,LocWt,LocAsize,LocPhi,[initialized]) variable LocDelta,LocNa,LocWt,LocAsize,Locphi variable initialized

//Default Parameters if(ParamisDefault(initialized)) initialized=0 endif

nvar omDelta,omNA,omf0,omASize,omPhi

omDelta=LocDelta; omNA=LocNa; omf0=LocWt; omASize=LocAsize; omPhi=LocPhi; if(! initialized) omRecalc("",0,"","") endif nvar omMass,omf0

variable f=sqrt((2*pi*omf0)^2+(avgWell(getKOpt))/omMass)/2/pi return f end //return to omFunctionList

//------------------------------------------------------------------------------------------------ // function getGQed() // Calculates the single-atom cavity QED coupling parameter in cyclic frequency // units. // // gQED = g0^2 / Dca //------------------------------------------------------------------------------------------------

function getGQed() nvar omNA, omG0, omDca return omG0^2/omDca end

//------------------------------------------------------------------------------------------------ // function getZHo() // Calculates the (single-atom) harmonic oscillator length // // zHO = sqrt(h / 2 m f0) //------------------------------------------------------------------------------------------------

function getZHo() nvar hbar, omMass, omf0 return sqrt(hbar/omMass/4/pi/omf0) end //return to omFunctionList

//------------------------------------------------------------------------------------------------ // function getLambDicke() // Calculates the Lamb-Dicke parameter // // LD = sqrt(2*nth+1)*zHO*kt //------------------------------------------------------------------------------------------------

function getLambDicke() nvar omNPhonon, omKT return sqrt(2*omNPhonon+1)*getZHo()*omKT end

//------------------------------------------------------------------------------------------------ // function getGranularity() // Calculates the (many-atom) granularity parameter using global // values of optomechanics parameters. // // Uses // epsilon = sqrt(N) F_dip z_HO / hbar kappa // where // F_dip = hbar kp g0^2 / dca // // c.f. Murch thesis //-------------------------------------------------------------------------------------------------

function getGranularity() nvar omKP, omK, omNA return omKP*getZHO()*sqrt(omNA)*getGQed() / omK end //return to omFunctionList

//------------------------------------------------------------------------------------------------- // function getHeating() // Returns the per-atom heating rate, in quanta / s // // Assumes heating only from Snn^(-) // // Heating is therefore 2*pi*nbar*kappa*epsilon^2/N //--------------------------------------------------------------------------------------------------

function getHeating(phi) variable phi nvar omK, omNBar, omNA variable ld ld = sqrt( sin(2*pi*phi)^2 + getLambDicke()^2*cos(2*pi*phi)^2 ) return 4*pi^2*omNBar*omK*getGranularity()^2/omNA*ld // omK is in cyclic units end //--------------------------------------------------------------------------------------------------- //function gOpt calculates OM damping rate function gOpt(freq) variable freq nvar omK, omf0

return 2 * omK * (omf0^2 - freq^2) / W1(freq) end //function W0 = kappa^2 + delta^2 function W0() nvar omK, omDelta

end

// function W1 = (kappa^2 + delta^2 - freq^2 ) function W1(freq) variable freq nvar omK, omDelta

return omK^2 + omDelta^2 - freq^2 end // function ksf not finished function ksf(freq) variable freq nvar omfs

return (omfs^2 - freq^2)* W0()/W1(freq) end

//--------------------------------------------------------------------------- //function FreqVsDelta_fit(w,t) : fitfunc // Fitting Function that calls GetOscFreq() // The x wave should be a list of deltas. // Uses cyclic frequencies for parameters //-----------------------------------------------------------------------------

function FreqVsDelta_fit(w,t) : fitfunc wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 6 //CurveFitDialog/ w[0] = omDelta_scaling //CurveFitDialog/ w[1] = omDelta_offset //CurveFitDialog/ w[2] = AtomNum //CurveFitDialog/ w[3] = Unshifted_TrapFreq //CurveFitDialog/ w[4] = AtomSize //CurveFitDialog/ w[5] = omPhi

return GetOscFreq(w[0]*t+w[1],w[2],w[3],w[4],w[5])

end //return to omFunctionList

//--------------------------------------------------------------------------- //function FreqVsPhi_fit(w,t) : fitfunc // Fitting Function that calls GetOscFreq() // The x wave should be a list of omPhi's // Uses cyclic frequencies for parameters //-----------------------------------------------------------------------------

function FreqVsPhi_fit(w,t) : fitfunc wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = omDelta //CurveFitDialog/ w[1] = AtomNum //CurveFitDialog/ w[2] = Unshifted_TrapFreq //CurveFitDialog/ w[3] = AtomSize //CurveFitDialog/ w[4] = omPhi_Offset

return GetOscFreq(w[0],w[1],w[2],w[3],t+w[4])

end

//--------------------------------------- // FITTING FUNCTIONS //---------------------------------------

//----------------------------------------------------------------------------------------------- //function omFitMaster(c,f,mode,[doLog]) // // Master function for all OM response fits. Has variable // functionality depending on choice of "mode" flag (see function // comments). For modes 0 and 1, returns a power spectral density // as a function of frequency "f". For mode 2 and 3, returns the real // and imaginary parts, respectively, of the response to transmitted // AM drive. // // Note: This function is NOT a fitting function! // Use a derived function. //----------------------------------------------------------------------------------------------- function omFitMaster(c,f,[mode,doLog]) wave c // c is the coefficients wave. Note that all frequencies must have identical units // (i.e. rad/s, kHz, etc.) // WARNING: In order to preserve the functionality of the various functions that // use omFitMaster, the order of these coefficients should NOT be changed // without changing all dependent functions. Consider adding new coefficients // at the END of the existing list. // c[0] -- Natural frequency // c[1] -- Shifted frequency // c[2] -- Natural damping rate (frequency units) // c[3] -- Cavity half-linewidth (frequency units) // c[4] -- Cavity detuning (frequency units) // c[5] -- Shot noise level (used only for modes 0 & 1) // c[6] -- Maximum technical noise level // For mode 0 -- Units of noise level // For modes 2-3 -- AM drive level // c[7] -- Minimum technical noise level, (noise level, used only for mode 0) // c[8] -- Noise peak quadrature angle (radians, used only for mode 0) // c[9] -- Detection quadrature angle (used only for modes 0 & 1) // c[10] -- Detection efficiency (used only for modes 0 & 1) variable f // f is the response frequency variable mode // Controls operation of omFitMaster // Possible values of mode (default is 0): // 0 - Normal operation, PSD driven by shot noise and constant technical noise // 1 - PSD driven by shot noise and technical noise in waves "techAm" and "techFm" // Due to subtleties of technical noise, this mode can not be used for quadratures other than 0 or pi/2. // 2 - Network analysis, real response, driven by AM noise before cavity // 3 - As above, imaginary response variable doLog // If true, results are returned in logarithm form, base 10. Default is 0.

if (ParamIsDefault(mode)) mode = 0; endif if (ParamIsDefault(doLog)) doLog = 0; endif

// Coefficient parameters variable f0 = c[0] // Natural resonance frequency variable fS = c[1] // Shifted resonance frequency variable gM = c[2] // Natural damping rate variable k = c[3] // Linewidth variable d = c[4] // Detuning variable aSN = c[5] // Shot noise power variable aTMax = c[6] // Maximum technical noise power variable aTMin = c[7] // Technical noise power offset variable qT = c[8] // Technical noise peak phase variable q = c[9] // Quadrature angle variable eDet = c[10] // Detection efficiency

// Derived parameters variable ksR = fS^2 - f0^2 // Real spring constant variable gOpt = 2*k*(f0^2 - f^2)/(k^2 + d^2 - f^2) // Optomechanical damping rate variable gTot = gM + gOpt // Total damping rate variable denomM = (fS^2 - f^2)^2 + f^2*gOpt^2 // Gain denominator magnitude

// Intracavity stuff variable GAmR = ((f0^2 - f^2)*(fS^2 - f^2) + f^2*gOpt^2) / denomM // Real part of the AM gain response to intracavity AM variable GAmI = (f0^2 - fS^2)*f*gOpt / denomM // Imaginary part of the AM gain response to intracavity AM variable GFmR = -(ksR/d)*(k*(fS^2 - f^2) + f^2*GOpt) / denomM // Real part of the FM gain response to intracavity AM variable GFmI = -(ksR/d)*(f*(fS^2 - f^2) + f*k*GOpt) / denomM // Imaginary part of the FM gain response to intracavity AM

// Shot noise stuff if (mode==0 || mode==1)

if (mode==1 && mod(q,pi/2)!=0) Abort "Error in omFitMaster: Mode 1 can only be used with AM or FM quadratures (you chose "+num2str(q)+" rad)." endif

// Extracavity stuff -- note that the following explicitly ignores the high-frequency filtering of the cavity Lorentzian!

variable iRot = 2*k^2/(k^2+d^2-f^2) // Rotation of drive to same phase quadrature (i.e. AM to AM) variable qRot = 2*k*d/(k^2+d^2-f^2) // Rotation of drive to out of phase quadrature (i.e. AM to FM)

variable MAmSnR = (GAmR*iRot - 1)*cos(q) + (GFmR*iRot - qRot)*sin(q) // Real response to AM shot noise variable MAmSnI = (GAmI*iRot)*cos(q) + (GFmI*iRot)*sin(q) // Imaginary response to AM shot noise variable MFmSnR = (GAmR*qRot)*cos(q) + (iRot + GFmR*qRot - 1)*sin(q) // Real response to FM shot noise variable MFmSnI = (GAmI*qRot)*cos(q) + (GFmI*qRot)*sin(q) // Imaginary response to FM shot noise

variable SSn = MAmSnR^2 + MAmSnI^2 + MFmSnR^2 + MFmSnI^2 // Total power response to shot noise

// Technical noise variable DAmT, DQT, SQT if (mode==0) DAmT = aTMin + aTMax*cos(qT) // Technical noise power in AM DQT = aTMin + aTMax*cos(qT-q) // Technical noise power at this quadrature SQT = max(0,DQT - DAmT) // Technical noise contribution from non-AM at this quadrature and frequency // Can not be less than 0 else // mode 1 wave techAm, techFm DAmT = techAm(f) if (mod(q,pi)==0) SQT = 0 else // Calculating FM quadrature SQT = techFm(f) endif endif variable SGTech = ( (GAmR*cos(q) + GFmR*sin(q))^2 + (GAmI*cos(q) + GFmI*sin(q))^2 )*DAmT // Atomic response to technical drive

variable STotal = (SSn + SGTech + SQT)*aSN // Total noise power STotal = eDet*STotal + aSN*(1-eDet) // Correct for detection efficiency

// Return output if (doLog) return log(STotal) else return STotal endif

endif

if (mode >= 2 ) // Network analysis stuff variable rotAmFm = (f^2*d^2)/((k^2 + d^2)^2 + f^2*k^2) variable SAm, SFm

if (mode==3) SAm = GAmR + 1 // Real AM to real AM SFm = GFmR + rotAmFm // Real AM to real FM elseif (mode==4) SAm = GAmI // Real AM to imaginary AM SFm = GFmI // Real AM to imaginary FM else abort "Error in omFitMaster. Mode was "+num2str(mode)+", must be between 0 and 3." endif

return (SAm*cos(q) + SFm*sin(q))*aTMax endif

end

//----------------------------------------------------------------------------------------------- // function makeOmPhaseMap(c,m) // // Makes a theoretical map of power response vs frequency and // quadrature phase. // // c is a coefficient wave containing: // c[0] -- Natural frequency // c[1] -- Shifted frequency // c[2] -- Natural damping rate (frequency units) // c[3] -- Cavity half-linewidth (frequency units) // c[4] -- Cavity detuning (frequency units) // c[5] -- Shot noise level // c[6] -- Maximum technical noise level (units of noise*frequency^2) // c[7] -- Minimum technical noise level (units of noise*frequency^2) // c[8] -- Noise peak quadrature angle (radians) // c[9] -- Detection efficiency // Here the technical noise is assumed to decay as 1/f^2 in all // quadratures. // // m is the destination wave. Frequency is indexed by rows, // and quadrature phase by columns. The axes must be properly // scaled before calling makeOmPhaseMap. //----------------------------------------------------------------------------------------------- function makeOmPhaseMap(c,m) wave c, m

// Make the coefficient wave for omFitMaster make /n=11 cTemp cTemp[0,8] = c[p] cTemp[10] = c[9]

variable j, i for (j=0; j<DimSize(m,1); j+=1) // Loop over quadrature angles cTemp[9] = DimOffset(m,1) + j*DimDelta(m,1) // This is the quadrature angle m[][j] = omFitMaster(cTemp,x,mode=0) endfor

killwaves/z cTemp end

//----------------------------------------------------------------------------------------------- // function omResponse(c,f) : fitfunc // // Fitting function for linear optomechanical power response // Coefficients are: // c[0] -- Natural frequency // c[1] -- Shifted frequency // c[2] -- Natural damping rate (frequency units) // c[3] -- Cavity half-linewidth (frequency units) // c[4] -- Cavity detuning (frequency units) // c[5] -- Shot noise level // c[6] -- Maximum technical noise level // c[7] -- Minimum technical noise level // c[8] -- Noise peak quadrature angle (radians) // c[9] -- Detection quadrature angle (radians) // c[10] -- Detection efficiency //----------------------------------------------------------------------------------------------- function omResponse(c,f) : fitfunc wave c variable f //CurveFitDialog/ //CurveFitDialog/ Coefficients 11 //CurveFitDialog/ c[0] -- f0 //CurveFitDialog/ c[1] -- fS //CurveFitDialog/ c[2] -- G //CurveFitDialog/ c[3] -- k //CurveFitDialog/ c[4] -- d //CurveFitDialog/ c[5] -- ASN //CurveFitDialog/ c[6] -- ATH //CurveFitDialog/ c[7] -- ATL //CurveFitDialog/ c[8] -- phiT //CurveFitDialog/ c[9] -- phi //CurveFitDialog/ c[10] -- eps return omFitMaster(c,f,mode=0) end

//----------------------------------------------------------------------------------------------- // function logOmResponse(c,f) : fitfunc // // Fitting function for linear optomechanical power response, // returning log of power // Coefficients are: // c[0] -- Natural frequency // c[1] -- Shifted frequency // c[2] -- Natural damping rate (frequency units) // c[3] -- Cavity half-linewidth (frequency units) // c[4] -- Cavity detuning (frequency units) // c[5] -- Maximum technical noise level (relative to shot noise) // c[6] -- Minimum technical noise level (relative to shot noise // c[7] -- Noise peak quadrature angle (radians) // c[8] -- Detection quadrature angle (radians) // c[9] -- Detection efficiency //----------------------------------------------------------------------------------------------- function logOmResponse(c,f) : fitfunc wave c variable f //CurveFitDialog/ //CurveFitDialog/ Coefficients 11 //CurveFitDialog/ c[0] -- f0 //CurveFitDialog/ c[1] -- fS //CurveFitDialog/ c[2] -- G //CurveFitDialog/ c[3] -- k //CurveFitDialog/ c[4] -- d //CurveFitDialog/ c[5] -- ATH //CurveFitDialog/ c[6] -- ATL //CurveFitDialog/ c[7] -- phiT //CurveFitDialog/ c[8] -- phi //CurveFitDialog/ c[9] -- eps

// Coefficient wave for omFitMaster as shot noise level = 1 make /free/n=11 c2 c2[0,4] = c[p] c2[5] = 1 c2[6,10] = c[p-1]

return log(omFitMaster(c2,f,mode=0)) end

//----------------------------------------------------------------------------------------------- //function Fano(w,t) : fitfunc // Fitting Function for shot-noise driven Fano resonance with // finite detection efficiency. //----------------------------------------------------------------------------------------------- function fano(w,t) : fitfunc wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = amp //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = ws //CurveFitDialog/ w[3] = w0 //CurveFitDialog/ w[4] = gamma

return w[0]*((1-w[1])+w[1]*((t^2-w[3]^2)^2+t^2*w[4]^2)/((t^2-w[2]^2)^2+t^2*w[4]^2))

end //return to omFunctionList

//----------------------------------------------------------------------------------------------- //function logFano(w,t) : fitfunc // Fitting Function for Log of Fano() // //----------------------------------------------------------------------------------------------- function logFano(w,t) : fitfunc wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = noise //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = ws //CurveFitDialog/ w[3] = w0 //CurveFitDialog/ w[4] = gamma

return Log((1-w[1])+w[1]*((t^2-w[3]^2)^2+t^2*w[4]^2)/((t^2-w[2]^2)^2+t^2*w[4]^2)+w[0])

end //return to omFunctionList

//----------------------------------------------------------------------------------------------- //function BreitWigner(w,t) : fitfunc // Fitting Function for Breit - Wigner resonance // //----------------------------------------------------------------------------------------------- function BreitWigner(w,t) : fitfunc wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 4 //CurveFitDialog/ w[0] = Amp //CurveFitDialog/ w[1] = y0 //CurveFitDialog/ w[2] = ws //CurveFitDialog/ w[3] = gamma

return w[1] + w[0]*w[2]^2*w[3]^2/( (t^2 - w[2]^2)^2 + t^2 * w[3]^2) End //----------------------------------------------------------------------------------------------- //function myLor(w,t) : fitfunc // Fitting Function for Lorenztain with independent amplitude and width // //----------------------------------------------------------------------------------------------- Function myLor(w,t) : FitFunc Wave w Variable t

//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog. //CurveFitDialog/ Equation: //CurveFitDialog/ f(t) = y0 + Amp* (gamma/2)^2 / ( (t - x0)^2 + (gamma/2)^2 ) //CurveFitDialog/ End of Equation //CurveFitDialog/ Independent Variables 1 //CurveFitDialog/ t //CurveFitDialog/ Coefficients 4 //CurveFitDialog/ w[0] = Amp //CurveFitDialog/ w[1] = y0 //CurveFitDialog/ w[2] = x0 //CurveFitDialog/ w[3] = gamma

return w[1] + w[0]* (w[3]/2)^2 / ( (t - w[2])^2 + (w[3]/2)^2 ) End

//function BreitWigner2(w,t) : fitfunc // Fitting Function for Breit - Wigner resonance with nonzero center frequency // //----------------------------------------------------------------------------------------------- function BreitWigner2(w,t) : fitfunc wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = Amp //CurveFitDialog/ w[1] = y0 //CurveFitDialog/ w[2] = ws //CurveFitDialog/ w[3] = gamma //CurveFitDialog/ w[4] = wOffset

return w[1] + w[0]*(w[2]-w[4])^2*w[3]^2/( ((t-w[4])^2 - (w[2]-w[4])^2)^2 + (t-w[4])^2 * w[3]^2) End

//function FourBreitWigner(w,t) : fitfunc // Fitting Function for Breit - Wigner resonance with nonzero center frequency // //----------------------------------------------------------------------------------------------- function FourBreitWigner(w,t) : fitfunc wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = Amp //CurveFitDialog/ w[1] = y0 //CurveFitDialog/ w[2] = ws //CurveFitDialog/ w[3] = gamma //CurveFitDialog/ w[4] = spacing

variable ws0 = w[2], ws1 = w[2]+w[4], ws2 = w[2] + 2*w[4], ws3 = w[2] - w[4]

variable bw0 = ws0^2*w[3]^2/( (t^2 - ws0^2)^2 + t^2 * w[3]^2) variable bw1 = ws1^2*w[3]^2/( (t^2 - ws1^2)^2 + t^2 * w[3]^2) variable bw2 = ws2^2*w[3]^2/( (t^2 - ws2^2)^2 + t^2 * w[3]^2) variable bw3 = ws3^2*w[3]^2/( (t^2 - ws3^2)^2 + t^2 * w[3]^2)

return w[1] + w[0]*(bw0 + 2*bw1 + 1.5*bw2 + 0.5*bw3)/5

End //-------------- //----------------------------------------------------------------------------------------------- //function BothBreitWigner(w,t) : fitfunc // Fitting Function for Breit - Wigner resonance // //----------------------------------------------------------------------------------------------- function BothBreitWigner(w,t) : fitfunc wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 6 //CurveFitDialog/ w[0] = AmpRed //CurveFitDialog/ w[1] = y0 //CurveFitDialog/ w[2] = ws //CurveFitDialog/ w[3] = gamma //CurveFitDialog/ w[4] = carrier //CurveFitDialog/ w[5] = AmpBlue

variable AmpRed = w[0] variable wsRed = w[2] variable carrier = w[4] variable wsBlue = 2*carrier - wsRed variable AmpBlue = w[5]

return w[1] + AmpRed*wsRed^2*w[3]^2/( (t^2 - wsRed^2)^2 + t^2 * w[3]^2) + AmpBlue*wsBlue^2*w[3]^2/( (t^2 - wsBlue^2)^2 + t^2 * w[3]^2) End //function BothFourBreitWigner(w,t) : fitfunc // Fitting Function for Breit - Wigner resonance with nonzero center frequency // //----------------------------------------------------------------------------------------------- function BothFourBreitWigner(w,t) : fitfunc wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 7 //CurveFitDialog/ w[0] = AmpRed //CurveFitDialog/ w[1] = y0 //CurveFitDialog/ w[2] = ws //CurveFitDialog/ w[3] = gamma //CurveFitDialog/ w[4] = spacing //CurveFitDialog/ w[5] = carrier //CurveFitDialog/ w[6] = AmpBlue

variable AmpRed = w[0] variable wsRed = w[2] variable carrier = w[5] variable wsBlue = 2*carrier - wsRed variable AmpBlue = w[6]

variable ws0Red = wsRed, ws0Blue = wsBlue variable ws1Red = wsRed+w[4], ws1Blue = wsBlue - w[4] variable ws2Red = wsRed + 2*w[4], ws2Blue = wsBlue - 2*w[4] variable ws3Red = wsRed - w[4], ws3Blue = wsBlue + w[4]

variable bw0Red = ws0Red^2*w[3]^2/( (t^2 - ws0Red^2)^2 + t^2 * w[3]^2) variable bw1Red = ws1Red^2*w[3]^2/( (t^2 - ws1Red^2)^2 + t^2 * w[3]^2) variable bw2Red = ws2Red^2*w[3]^2/( (t^2 - ws2Red^2)^2 + t^2 * w[3]^2) variable bw3Red = ws3Red^2*w[3]^2/( (t^2 - ws3Red^2)^2 + t^2 * w[3]^2)

variable bw0Blue = ws0Blue^2*w[3]^2/( (t^2 - ws0Blue^2)^2 + t^2 * w[3]^2) variable bw1Blue = ws1Blue^2*w[3]^2/( (t^2 - ws1Blue^2)^2 + t^2 * w[3]^2) variable bw2Blue = ws2Blue^2*w[3]^2/( (t^2 - ws2Blue^2)^2 + t^2 * w[3]^2) variable bw3Blue = ws3Blue^2*w[3]^2/( (t^2 - ws3Blue^2)^2 + t^2 * w[3]^2)

variable amp0 = 1, amp1 = 2, amp2 = 1.4, amp3 = 0.6; variable ampSum = amp0 + amp1 + amp2 + amp3

return w[1] + AmpRed*(amp0*bw0Red + amp1*bw1Red + amp2*bw2Red + amp3*bw3Red)/ampSum + AmpBlue*(amp0*bw0Blue + amp1*bw1Blue + amp2*bw2Blue + amp3*bw3Blue)/ampSum

End //------------- //--------------------------------------------------------------------------------------------------- function totalAmResponse5(inx) variable inx wave w = w_coef nvar omTfreq variable t = omTfreq //CurveFitDialog/ //CurveFitDialog/ Coefficients 7 //CurveFitDialog/ w[0] = snLevel //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = fs //CurveFitDialog/ w[3] = f0 //CurveFitDialog/ w[4] = gammaM //CurveFitDialog/ w[5] = kappa //CurveFitDialog/ w[6] = Delta variable deltanoise = w[7]

variable ReGNum, ImGNum, GDenom, ReG, ImG, A, B, KD, k, d, fs, f0, gM variable techPwr wave techwaveAM

fs = w[2]; f0 = w[3]; k = w[5]; d = w[6]+inx; gM = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt

ReGNum = (fs^2 - f0^2)*(t^2 - fs^2) // Numerator of real part of gain function ImGNum = t*gTot*(fs^2-f0^2) // Numerator of imaginary part of gain function GDenom = (t^2 - fs^2)^2 + t^2*gTot^2 // Denominator (real) of gain function

ReG = ReGNum/GDenom ImG = ImGNum/GDenom

A = 2*k^2/(k^2 + d^2) // Ratio of 2*kappa^2 to kappa^2 + Delta^2 B = (k^2-d^2)/(k^2+d^2) // Unperturbed reflected field

// Real part of the shot noise response // Corrected by NB Feb 28, 2011 variable ReSn = ( ReG*A + B )^2 + ( ImG*A*d/k )^2 variable ImSn = ( ImG*A )^2 + ( (ReG+1)*A*d/k )^2 variable Sam if(w[0]==0) Sam = (ReG+1)^2 + ImG^2 else techPwr = (techWaveAM(t)/w[0] - 1)/w[1]

variable Tech = techPwr*( (ReG+1)^2 + ImG^2 )

Sam= w[0]*( w[1]*(ReSn + ImSn + Tech) + 1 - w[1] ) // Total AM optomechanical response endif return 1/sqrt(2*pi)/(DeltaNoise)*Exp(-(inx^2)/(2*DeltaNoise^2))*Sam End //return to omFunctionList

//return to omFunctionList //----------------------------------------------------------------------------------------------- //function totalFmResponse4(w,t) : fitfunc // Fitting Function for shot-noise driven BW resonance with // finite detection efficiency and optomechanical damping. // Corrected for phase error in Gain //----------------------------------------------------------------------------------------------- function totalFmResponse4(w,t) wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 7 //CurveFitDialog/ w[0] = snLevel //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = fs //CurveFitDialog/ w[3] = f0 //CurveFitDialog/ w[4] = gammaM //CurveFitDialog/ w[5] = kappa //CurveFitDialog/ w[6] = Delta

variable ReGNum, ImGNum, GDenom, ReG, ImG, A, B, KD, k, d, fs, f0, gM variable techPwr wave techwaveAM wave techwaveFM

fs = w[2]; f0 = w[3]; k = w[5]; d = w[6]; gM = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt

ReGNum = (fs^2 - f0^2)*(k*(fs^2 - t^2) + t^2*gtot) // Numerator of real part of gain function ImGNum = t*(fs^2-f0^2)*(k*gtot - (fs^2-t^2)) // Numerator of imaginary part of gain function GDenom = d*(fs^2 - t^2)^2 + d*t^2*gTot^2 // Denominator (real) of gain function

ReG = ReGNum/GDenom ImG = ImGNum/GDenom

A = 2*k^2/(k^2 + d^2) // Ratio of 2*kappa^2 to kappa^2 + Delta^2

variable ReSn = 2*A*(ReG^2 + ImG^2 + d/k*ReG)

variable aToA = ((k^2 + d^2)^2 + t^2*k^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) // AM to AM transduction variable aToF = (t^2*d^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) variable ratioFA = aToF/aToA

if(w[0]==0) return ReG^2 + ImG^2 + ratioFA else techPwr = (techWaveAM(t)/w[0] - 1)/w[1] variable Tech = techPwr*(ReG^2 + ImG^2)

return w[0]*( w[1]*(ReSn+1+Tech) + 1 - w[1] ) + techwaveFM(t) -w[0] // Total FM optomechanical response endif End //return to omFunctionList function fitavgFmResp5(w,t) : fitfunc wave w variable t

//CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = snLevel //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = fs //CurveFitDialog/ w[3] = f0 //CurveFitDialog/ w[4] = gammaM //CurveFitDialog/ w[5] = kappa //CurveFitDialog/ w[6] = Delta //CurveFitDialog/ w[7] = DeltaNoise

duplicate/o w wtemp wtemp=w variable deltaNoise = w[7] variable /g omTfreq = t return Integrate1d(totalFmResponse5,-4*deltaNoise,4*deltaNoise,0,100)

end

function totalFmResponse5(inx) variable inx wave w = w_coef nvar omTfreq variable t = omTfreq //CurveFitDialog/ //CurveFitDialog/ Coefficients 7 //CurveFitDialog/ w[0] = snLevel //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = fs //CurveFitDialog/ w[3] = f0 //CurveFitDialog/ w[4] = gammaM //CurveFitDialog/ w[5] = kappa //CurveFitDialog/ w[6] = Delta variable deltaNoise = w[7]

variable ReGNum, ImGNum, GDenom, ReG, ImG, A, B, KD, k, d, fs, f0, gM variable techPwr wave techwaveAM wave techwaveFM

fs = w[2]; f0 = w[3]; k = w[5]; d = w[6] + inx; gM = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt

ReGNum = (fs^2 - f0^2)*(k*(fs^2 - t^2) + t^2*gtot) // Numerator of real part of gain function ImGNum = t*(fs^2-f0^2)*(k*gtot - (fs^2-t^2)) // Numerator of imaginary part of gain function GDenom = d*(fs^2 - t^2)^2 + d*t^2*gTot^2 // Denominator (real) of gain function

ReG = ReGNum/GDenom ImG = ImGNum/GDenom

A = 2*k^2/(k^2 + d^2) // Ratio of 2*kappa^2 to kappa^2 + Delta^2

variable aToA = ((k^2 + d^2)^2 + t^2*k^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) // AM to AM transduction variable aToF = (t^2*d^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) variable ratioFA = aToF/aToA

variable ReSn = 2*A*(ReG^2 + ImG^2 + d/k*ReG) variable Sfm

if(w[0]==0) Sfm = ReG^2 + ImG^2 + ratioFA else techPwr = (techWaveAM(t)/w[0] - 1)/w[1] variable Tech = techPwr*(ReG^2 + ImG^2)

Sfm = w[0]*( w[1]*(ReSn+1+Tech) + 1 - w[1] ) + techwaveFM(t) -w[0] // Total FM optomechanical response endif

return 1/sqrt(2*pi)/(DeltaNoise)*Exp(-(inx^2)/(2*DeltaNoise^2))*Sfm End //return to omFunctionList //----------------------------------------------------------------------------------------------- //function AMresponseAMdrive(w,t) : fitfunc // Fitting Function for AM technical noise driven data // Plots the response in the AM quadrature of // (Amplitude_Atoms - Amplitude_SN)/(Amplitude_Empty - Amplitude_SN) //----------------------------------------------------------------------------------------------- function AMresponseAMdrive(w,t) wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs //CurveFitDialog/ w[1] = f0 //CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta

Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt

Variable ResponseNum = ((fs^2-t^2)*(f0^2-t^2)+(t*gTot)^2)^2+(t*gTot*(fs^2-f0^2))^2 Variable ResponseDenom = ((fs^2-t^2)^2+(t*gTot)^2)^2

Return ResponseNum/ResponseDenom

End //return to omFunctionList

//----------------------------------------------------------------------------------------------- //function FMresponseAMdrive(w,t) : fitfunc // Fitting Function for AM technical noise driven data // Plots the response in the FM quadrature of // (Amplitude_Atoms - Amplitude_SN)/(Amplitude_Empty - Amplitude_SN) //----------------------------------------------------------------------------------------------- function FMresponseAMdrive(w,t) wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs //CurveFitDialog/ w[1] = f0 //CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta

Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt

Variable RealNum = d^2*t*((fs^2-t^2)+(t*gTot)^2)+k*(d^2+k^2)*t*gTot*(fs^2-f0^2) Variable ImagNum = -k*(d^2+k^2)*(fs^2-f0^2)*(fs^2-t^2) Variable ResponseDenom = d*(d^2+k^2)*((fs^2-t^2)^2+(t*gTot)^2)

variable aToA = ((k^2 + d^2)^2 + t^2*k^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) // AM to AM transduction variable aToF = (t^2*d^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) variable ratioFA = aToF/aToA

Return (RealNum^2+ImagNum^2)/(ResponseDenom)^2 + ratioFA

End //return to omFunctionList

//----------------------------------------------------------------------------------------------- //function AMIPresponseAMdrive(w,t) : fitfunc // Fitting Function for the normalized in-phase response in the // AM quadrature to an AM drive // The function fits to the following curve // (Note: Amp = Amplitude and SN = Shot-noise) // (Amp_Atoms-Amp_SN)/(Amp_EmptyCavity-Amp_SN) //----------------------------------------------------------------------------------------------- function AMIPresponseAMdrive(w,t) wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs //CurveFitDialog/ w[1] = f0 //CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta

Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM+gOpt

Variable ReG = ((f0^2-fs^2)*(fs^2-t^2))/((fs^2-t^2)^2+(t*gTot)^2)

Return 1+ReG

End //return to omFunctionList

//----------------------------------------------------------------------------------------------- //function AMOPresponseAMdrive(w,t) : fitfunc // Fitting Function for the normalized out-of-phase response in the // AM quadrature to an AM drive // The function fits to the following curve // (Note: Amp = Amplitude and SN = Shot-noise) // Log((Amp_Atoms-Amp_SN)/(Amp_EmptyCavity-Amp_SN)) //----------------------------------------------------------------------------------------------- function AMOPresponseAMdrive(w,t) wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs //CurveFitDialog/ w[1] = f0 //CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta

Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt

Variable ImG = (f0^2-fs^2)*(t*gTot)/((fs^2-t^2)^2+(t*gTot)^2)

Return -ImG

End //return to omFunctionList

//----------------------------------------------------------------------------------------------- //function FMresponseAMdrive(w,t) : fitfunc // Fitting Function for the normalized in-phase response in the // FM quadrature to an AM drive // The function fits to the following curve // (Note: Amp = Amplitude and SN = Shot-noise) // (Amp_Atoms-Amp_SN)/(Amp_EmptyCavityAM-Amp_SN) //----------------------------------------------------------------------------------------------- function FMIpResponseAmDrive(w,t) wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs //CurveFitDialog/ w[1] = f0 //CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta

Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt

Variable ReG = ((f0^2-fs^2)*(fs^2-t^2))/((fs^2-t^2)^2+(t*gTot)^2)

variable ratioReFA = d*t/sqrt((k^2+d^2)^2 + t^2*k^2) // May need a - sign prefactor depending on Fourier convention

Return -k*ReG/d + ratioReFA

End //return to omFunctionList

//----------------------------------------------------------------------------------------------- //function FMresponseAMdrive(w,t) : fitfunc // Fitting Function for the normalized out-of-phase response in the // FM quadrature to an AM drive // The function fits to the following curve // (Note: Amp = Amplitude and SN = Shot-noise) // Log((Amp_Atoms-Amp_SN)/(Amp_EmptyCavityAM-Amp_SN)) //----------------------------------------------------------------------------------------------- function FmOpResponseAMdrive(w,t) wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs //CurveFitDialog/ w[1] = f0 //CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta

Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt

// This is what Thierry calculates Variable ImG = (f0^2-fs^2)*(t*gTot)/((fs^2-t^2)^2+(t*gTot)^2)

// This is what I think it may be, but is from what I used to calculate the SN response // this has extra term from -i*w in k - i*w in numerator:

// Variable ImG = t*(fs^2 - f0^2)*(gTot - (fs^2 - t^2)/k) / ( (fs^2- t^2)^2 + t^2*gTot^2)

Return k*ImG/d

End //return to omFunctionList

function AllphaseAllquadResponseAmDrive(w,t) wave w variable t nvar npoints //CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = fs //CurveFitDialog/ w[1] = f0 //CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta

wave rAmIp variable deltaFreq = npoints*deltax(rAmIp)

if(t< (leftx(rAmIp) + deltaFreq)) return AMIPresponseAMdrive(w,t) elseif(t< (leftx(rAmIp) + 2*deltaFreq)) return AMOPresponseAMdrive(w,t-deltaFreq) elseif(t< (leftx(rAmIp) + 3*deltaFreq)) return FMIPresponseAMdrive(w,t-2*deltaFreq) else return FMOPresponseAMdrive(w,t-3*deltaFreq) endif end

function makeRelief34(zmin,zmax) variable zmin, zmax

colortab2wave relief

duplicate/o m_colors ctab setscale /i x,zmin,zmax,ctab

variable i for (i=126; i>=0; i-=2) deletepoints (i),1,ctab endfor

end

//return to omFunctionList

//----------------------------------------------------------------------------------------------- //function SSresponse (w,t) : fitfunc // Fitting Function for the single-sideband on-resonance om // response // //----------------------------------------------------------------------------------------------- function SSresponse(w,t) : fitfunc wave w variable t //CurveFitDialog/ //CurveFitDialog/ Coefficients 10 //CurveFitDialog/ w[0] = gom //CurveFitDialog/ w[1] = wm //CurveFitDialog/ w[2] = Gamma_m //CurveFitDialog/ w[3] = p_th //CurveFitDialog/ w[4] = SN_Amplitude //CurveFitDialog/ w[5] = Offset //CurveFitDialog/ w[6] = Het_Det_Efficiency //CurveFitDialog/ w[7] = Kappa //CurveFitDialog/ w[8] = w_cavity //CurveFitDialog/ w[9] = nbar

variable gom = w[0] variable wm = w[1] variable Gm = w[2] variable pth = w[3] variable SNamp = w[4] // Ideally, this should equal P_LO variable SNoffset = w[5] // Ideally, this should equal 1 variable EffHet = w[6] variable kappa = w[7] variable fc = w[8] // frequency of cavity resonance variable nbar = w[9]

variable Copt = 4*nbar*gom^2 / kappa / Gm

variable denom, numerator denom = (wm^2+(Gm/2)^2-(t-fc)^2)^2 + ((t-fc)*Gm)^2 numerator = 2 * wm * (Copt*wm-(t-fc)) + (wm^2 + (Gm/2)^2+(t-fc)^2) * (2*pth+1) return SNamp *(SNoffset + EffHet*Copt*Gm^2 * numerator / denom) End

//return to omFunctionList