Some clarification is warranted regarding the different shapes of stacks. Corrections are appreciated:

1) The single-leg stack that I just completed should function as a working model for the IO, OO, and MC1/3. Rana commented, however, that the current dimensions are slightly off for MC1/3 (which makes sense since I could only find drawings for the IOC). If anyone knows the whereabouts of similar drawings for MC1/3, I'd much appreciate it.

2) A triple-leg stack can model the BS, ITMX, and ITMY chambers. I don't have exact dimensions for these, but I can make decent approximations from to-scale stack drawings. I'll probably work on a model for this style next, since at least I have some information regarding this version.

3) The MC2 chamber has its own stack model, about which I haven't found any drawings in the binders. I can't start on MC2C at all until I find such drawings.

7/14: Analytical calculation of Viton spring constant; updated Viton values in models; experimental confirmation of COMSOL eigenfrequencies (single stack layer)

7/15: Extensions to 2-, 3-, and 4-layer stack legs. Eigenfrequency characterizations performed for each level. Meshing issues with 4-layer stack prevented completion.

7/19: Debugged the 4-layer stack. Turned out to be a boundary condition issue because of non-sequential work-plane definitions. Successful characterization of single-leg eigenfrequencies.

7/20: Prototype three-legged stack completed, but dimensions are incorrect. Read Sievers paper for details of triple-legged stack. Sorted through many stack design binders in efforts to distinguish IOC/OOC, BSC/ITMX/ITMY, MC1/MC3, and MC2 dimensions.

7/21: Researched frequency domain analysis testing in COMSOL. Attempting to first find transfer function of a single-layer stack --> currently running into some run-time errors that will need some more debugging in the afternoon.

Now that venting is complete, this is a request for anyone who opens any chamber:

1) Please notify me immediately so I can take pictures of the stacks in that chamber.

2) If I'm not around, please take a few pictures for me. I'm most interested in the shape, number of layers, size, and damper arrangements of each stack.

This is most important for the MC1/MC3 chamber, MC2 chamber, and BS/ITMX/ITMY chambers.

Over the past couple days, I discovered a simple, direct method for calculating frequency responses with a combination of COMSOL and any plotter such as Excel or MatLab. The simple case of rectangular prism of steel was analyzed using this method; details will be posted shortly on the COMSOL Wiki page. The frequency response matched theoretical reasoning: the bar acts as a simple mechanical low-pass filter, rapidly attenuating driving frequencies at the base beyond the first eigenmode.

It therefore shouldn't be too difficult to extend this analysis to the MC1/MC3 stack. The many eigenfrequencies will produce a more complicated transfer function, and so more data points will be taken.

The major shortcoming of this method involves dealing with the imaginary components of the eigenfrequencies. As of now, I haven't found a way of measuring the phase lag between the drive and the response. I also haven't found a way of changing the damping constants and therefore playing with phase components.

I have successfully completed a preliminary transfer function measurement test on the MC1/MC3 stack in COMSOL. Using the measurement scheme described on the Wiki, I initialized a 1 N/m^2 sinusoidal perturbation on the bottom of the stack and measured the maximum displacement of the top layer. This preliminary test just calculated the responses to 1-,2-,3-,4-, and 5-Hz drives along the x-axis (pictures attached).

Currently, I am rerunning the same test but from 1-10 Hz with 0.1-Hz steps. When both x- and y-axis responses have been plotted, I can move on to repeating this entire process on the MC2 stack.

I completed the frequency domain analysis mentioned previously in the x-direction. Although I ran it from 1-10 Hz, with 0.1-Hz increments, COMSOL was unable to complete the task past 7 Hz because the relative error was beyond the relative tolerance. To solve this issue, I'd have to rerun the simulation with a finer mesh, an unfavorable option because of the already-extensive run times. The Bode magnitude plot from this simulation is attached:

This simulation raises some questions about the feasibility of this method:

1) Do we have the computing power necessary?

I already moved my work from my personal Mac Pro to Kallo (4 GB --> 12 GB RAM difference). Now, instead of crashing the program constantly, I typically wait a half hour for a standard run of the model. Preferably, I could move my work to Megatron or some other workhorse-computer... but I also know that many of the big boys are already being strained as is.

2) Is it possible to take measurements through Matlab?

This way, I could write a script to instruct COMSOL and just run a few tests at a time overnight. Also, I wouldn't have to sit and record measurements manually as I've done here. The benefits of such an improvement warrant further attention. I'll work on this option next.

3) Up until what frequency do we need to model?

As I've shown, normal meshing yields data up to 7 Hz. Is this enough? Do we need more data? Certainly not less, I'm quite sure. We need to keep in mind that as frequency range increases, run times increase exponentially.

4) How do we incorporate gravity into the equation?

Gravity will produce a bit of extra force in the non-restoring direction for off-axis deviations, slightly decreasing the expected frequency. Whether or not this is an important effect is questionable, since the deviations are typically on the order of a micron, which is orders of magnitude smaller than the initial displacement I'm using on the base. I've decided to ignore this complication for now.

1) Gravity has to be included because the inverted pendulum effect changes the resonant frequencies. The deflection from gravity is tiny but the change in the dynamics is not. The results are not accurate without it. The z-direction probably is unaffected by gravity, but the tilt modes really feel it.

2) You should try a better meshing. Right now COMSOL is calculating a lot of strain/stress in the steel plates. For our purposes, we can imagine that the steel is infinitely stiff. There are options in COMSOL to change the meshing density in the different materials - as we can see from your previous plots, all the action is in the rubber.

3) I don't think the mesh density directly limits the upper measurement frequency. When you redo the swept-sine using the matlab scripting, use a logarithmic frequency grid like we usually do for the Bode plots. The measurement axis should go from 0.1 - 30 Hz and have ~100 points.

In any case, the whole thing looks promising: we've got real solid models and we're on the merge of being able to duplicate numerically the Dugolini-Vass-Weinstein measurements.

I made some progress on a couple issues:

1) I figured out how to create log-transfer function plots directly in COMSOL, which eliminates the hassle of toggling between programs.

2) Instead of plotting maximum displacement, which could lead to inconsistencies, I've started using point displacement, standardizing to the center of the top surface.

3) I discovered that the displacement can be measured as a field vector, so the minor couplings between each translational direction (due to the asymmetry in the original designs) can be easily ignored.

For the past couple of days, Jan and I have been discussing a major issue in COMSOL involving modeling both oscillatory and non-oscillatory forces simultaneously while using FDA. It turns out that he and I had run into the same problem at different times and with different projects. After discussing with an expert, Jan had decided in the past that this simple task was impossible via direct means.

The issue could still be resolved if there was a way for us to work on the Weak Form of the differential equations describing the system:

Usually, one must define weight as a body load in the negative-z direction. However, this problematically instantiates a new force in COMSOL, which is automatically driven over the range of frequencies during FDA.

Instead, we could define gravity as an anti-restoring force, since we assume that the base of the stack is fixed.

In other words,F_{g} = (ρ*g/L)*x + (ρ*g/L)*y for a point mass which is constrained on the bottom (for small angles).

Working in Weak Form then, we'd never have to define an explicit gravity load-- this could just be an extra couple of terms in the differential equation which are related entirely to the x- and y-vectors (well-defined for each mesh point). This would fool COMSOL into never tacking on the oscillatory term during FDA.

According to current documentation however, Weak Form analysis is not yet possible in COMSOL 4.0. Jan suggested moving my work over to ANSYS or waiting for the 4.0 upgrade, but there's probably not enough time left in my SURF for either of these options. I suggested attempting a backwards-compatibility test to COMSOL 3.5; Jan and I will be exploring this option some time next week.

I have discovered a method of completely characterizing the 6x6 response of all six types (x-,y-, and z- translational/rotational) of oscillatory disturbances at the base of the stack.

"Tipping" drives are trivial, and simply require a face load in the appropriate direction.

"Tilting" drives could use a torque, but I am instead implementing multiple edge loads in opposing directions to create the appropriate net curl. This curl will be kept constant across the three axes for sake of comparing the resulting transfer functions.

"Tipping" responses are once again trivial, and merely require the displacement vector of the top center coordinate to be recorded.

"Tilting" responses require the normal vector to be recorded and manipulated to produce the angular coordinates (assuming right-handed coordinate system):

θ_{x }= tan^{-1}(x/z)

θ_{y }= tan^{-1}(y/z)

θ_{z }= tan^{-1}(y/x)

The first three concepts have been confirmed through simulations to produce correct transfer functions. The last test seems to be producing some problems, in that the vector normal to the equilibrium position (an obvious and useless piece of information) is sometimes given instead of the vector normal to the position of maximum displacement. This means that, as of now, I have the capability of measure the half of the complete 6x6 matrix of transfer functions in the coming weeks. The first three of eighteen transfer functions are attached below and will be included in my progress report.

I reran the FDA in COMSOL on the MC1/MC3 Stack and produced the following Displacement-Displacement Transfer Functions:

X-Translational Drive has a blue background

Y-Translational Drive has a red background

Z-Translational Drive has a green background

Obtaining the Displacement-to-Phase part of the Transfer Function still produces difficulties -- I'm still working on the COMSOL-Matlab interface to perhaps better facilitate this.

RA: I have deleted those plots because they weren't transfer functions. Transfer functions must always be the ratio of something to something. For example: if I had a nickel for every bad plot I see, I would be a millionaire. In that example, the transfer function would have the units of nickels/plots. For the stacks, it should be meters/meter.

I reran the FDA in COMSOL on the MC1/MC3 Stack and produced the following Displacement-Displacement Transfer Functions:

X-Translational Drive has a blue background

Y-Translational Drive has a red background

Z-Translational Drive has a green background

Obtaining the Displacement-to-Phase part of the Transfer Function still produces difficulties -- I'm still working on the COMSOL-Matlab interface to perhaps better facilitate this.

RA: I have deleted those plots because they weren't transfer functions. Transfer functions must always be the ratio of something to something. For example: if I had a nickel for every bad plot I see, I would be a millionaire. In that example, the transfer function would have the units of nickels/plots. For the stacks, it should be meters/meter.

My apologies for the mislabeled axes on my previous plots. They have been corrected to a ratio (in./in.), as Rana so kindly suggested in his helpful, not-at-all-condescending response.

I have chosen to stay in the English system because all of the original stack drawings are in inches as well.

Time Domain Analysis on a Driven, Damped Simple Pendulum has produced a method for implementing gravity.

COMSOL made this simple task a cryptic one: the following methods had previously failed:

Previous Frequency Domain testing lead to unwanted oscillations of all loads.

Prescribed accelerations at first seemed to create a constant gravity, but instead incorrectly constrained net acceleration to the inputted amount

Methodology:

1) An (approximately) impulse displacement was applied in the horizontal direction. The pendulum bob's displacement was observed for varying pendulum lengths.

2) The drive and response displacements vs. time were FFT'd to produce transfer functions.

3) The fundamental frequencies were inverted, squared, and plotted against frequency.

4) Since the graph is linear with an R^2 of over 0.99, it is reasonable to assume that gravity is properly acting as a restoration force.

For the sake of future users, I have decided to periodically add tips and tricks in using COMSOL that I have figured out, most probably after hours of circuitous efforts. They will always be listed under the new COMSOL Tips category.

Today's topic: Intrusions

COMSOL has a very user-friendly interface for taking objects from 2D to 3D using the "extrusion" feature. But suppose one wants to design an object which contains screw holes or some other indentation. I've found that creating "punctures" in COMSOL is either impossible or very complicated.

Instead, COMSOL encourages users to always "add" to the object. In other words, one must form the lowest level first, then build layers sequentially on top using new work plane and boolean difference operators. This will probably be a bit clearer with an example:

1) First, create the planar projection in a work plane:

2) Extrude the first layer only in the regular fashion:

3) Add a new work plane which is offset in the z-direction to the deepest point of the intrusion.

4) Now, create the shape of the intrusion in this new work plane.

5) Use the Boolean "Difference" to let COMSOL know that, on this plane, the object has a hole.

6) Extrude once more from the second work plane to complete the intrusion.

Wed. 7/7: COMSOL Busbar tutorials; began stack design; began base; Viton rubber research

Thurs. 7/8: Completed Viton rubber research; updated materials; finished designing the base layer

Fri. 7/9: Research model coupling papers; extensive eLog entry about base design and troubleshooting

Sun. 7/11: Played around with Busbar to find first eigenfrequency; continued crashing COMSOL

Mon. 7/12: Intrusions in COMSOL eLog tutorial entry; research eigenfrequency analysis; successfully got first eigenmode of rectangular bar

Tues. 7/13: Updated Poisson ratio of Viton and subsequently succeeded in running eigenfrequency tests on base stack layer. Systematic Perturbation Tests were documented in the most recent elog entry. Discussed results with Rana and decided this didn't make sense. Analytical study required.

Wed. 7/14: Went over to machine shop to experimentally extrapolate spring constant of Viton. Calculations to be done in the afternoon.

In this experiment, we used a feedback control to create a stable trap for a NdFeB permanent magnet. The block diagram is the following:

The displacement of the magnet is sensed by the Hall-effect sensor, whose output voltage is proportional to the magnetic flux produced

by the permanent magnet. It has a flat response within the frequencies we are interested in. It is driven by a 5 V power supplier and its

output has a DC voltagle of 2.5 V. We subtracted the DC voltage and used the resulting signal as the error signal. This was

simply achieved by using two channels "A" and "B". The output is "A-B" with a gain equal to one. We then put the error

signal into another SR560 as a low-pass filter with a gain of 100 above 30 Hz. We used the "DC" coupling modes in both

preamplifers. The output is then used to drive a coil. The coil has a dimension of 1.5 inch in diameter and 2 inch in length.

The inductance of the coil is around 0.5 H and the resistance is 4.7 Om. Therefore, it has a corner frequency aournd 10/2pi Hz.

The coil has a iron core inside to enhance the DC force to the permanent magnet. The low-frequency 1/f response of the magnet is produced

by the eddy current damping of the aluminum plane that is below the magnet. This 1/f response is essential for a stable configuration. In the

next stage, we will remove the aluminum plane, and instead we will use a filter to create similar transfer function. At high-frequencies, it behaves as

a free-mass and has a 1/f^2 response. Finally, the magnet is stably levitated.

I am using SR785 Spectrum Analyzer now and also tomorrow.
I will put it back on Sunday. If anyone needs it during the weekend,
please just take it and I can reset it by myself later. Thanks.

I measured the open-loop transfer function of the magnetic levitation system.

The schematic block diagram for this measurement is the following:

I injected a signal at a level of 20mV between two preamplifiers, and the corresponding open-loop

transfer function is given by B/A. I took a picture of the resulting measurement, because

I encountered some difficulties to save the data to the computer via the wireless network.

The bode plots for the transfer function shown on the screen is the following:

I am puzzled with the zero near 10 Hz. I think it should come from the mechanical response function, because there is no zero in the transfer functions

of the preamplifer and the coil itself. I am not sure at the moment.

The corresponding configuration of the levitated magnet is

1. Why do all of the BNCs have no GND connection? Each should have the individual cables to the ground. Each signal line and the corresponding ground line should be twisted together.

2. This looks the (usual) oscillation of the V-I conversion stage but I can't tell anything as I don't have the circuit diagram of the whole circuit.

3. In a certain case, putting some capacitance at the feedback may help. Read P.11 of the data sheet of LT1125. Try to put some capacitors from 20pF to some larger one whether it changes the situation or not. I suppose the bandwidth of your sensor can be ~1kHz. So putting a capacitance less than 10nF still has no effect to the servo.

1. They are all connected to the box which has a single connection to the board ground. If I connect each of them to the ground, there would be many small loops

of ground. Of course, I should have replaced all the connectors such that the they are disconnected to the box as point out by Robert.

2. The oscillation disappears after I add 5 nF capacitor in parallel to the existing resistor. Thank you very much for pointing this out.

We did some quick DC balancing of the MC2 coil drivers to reduce the l2a coupling. We updated the gains in the C1:SUS-MC2_UL/UR/LR/LLCOIL to be 1, -0.99, 0.937,-0.933, respectively. The previous values were 1, -1, 1, -1.

The procedures are the following:

Lock IMC.

Drive UL+LR and change the gain of LR to zero pitch.

Drive UR+LL and change the gain of LL to zero pitch.

Lastly, drive all 4 coils and change UR & LR together to zero yaw.

We used C1:SUS-MC2_LOCKIN1_OSC to create the excitations at 33 Hz w/ 30,000 cts. The angular error signals were derived from IMC WFSs.

While this time we did things by hand, in the future it should be automated as the procedure is sufficiently straightforward.

My old scheme was flawed as I used pitch as the readback. The pitch signal could not distinguish the cross-coupling due to coil imbalance and that due to the natural suspension L2P. A new scheme based on yaw alone has been developed and will be integrated into ifo_test. For now we revert the C1:SUS-MC2_UL/UR/LR/LLCOIL gains back to 1, -1, 1, -1.

Quote:

[Yehonathan, Hang]

We did some quick DC balancing of the MC2 coil drivers to reduce the l2a coupling. We updated the gains in the C1:SUS-MC2_UL/UR/LR/LLCOIL to be 1, -0.99, 0.937,-0.933, respectively. The previous values were 1, -1, 1, -1.

The procedures are the following:

Lock IMC.

Drive UL+LR and change the gain of LR to zero pitch.

Drive UR+LL and change the gain of LL to zero pitch.

Lastly, drive all 4 coils and change UR & LR together to zero yaw.

We used C1:SUS-MC2_LOCKIN1_OSC to create the excitations at 33 Hz w/ 30,000 cts. The angular error signals were derived from IMC WFSs.

While this time we did things by hand, in the future it should be automated as the procedure is sufficiently straightforward.

After updating the 40 m finesse file to incorporate the new SRC length (and the removal of SR2), we find that the current SRM radius curvature is fine. Thus a replacement of SRM is NOT required.

Basically, the new one-way SRC gouy phase is 11.1 deg according to Finesse, which is very close to the previous value of 10.8 deg. Thus the transmode spacing should be essentially the same.

In the first attached plot is the mode content calculated with Finesse. Here we have first offset DARM by 1m deg and misaligned the SRM by 10 urad. From the top to bottom we show the amplitude of the carrier fields, f1, and f2 sidebands, respectively. The red vertical line is the nominal operating point (thanks Koji for pointing out that we do signal recycling instead of extraction now). No direct co-resonance for the low-order TEM modes. (Note that the HOMs appeared to also have peaks at \phi_srm = 0. This is just because the 00 mode is resonant and thus the seed for the HOMs is greater. )

We can also consider a clean case without mode interactions in the second plot. Indeed we don't see co-resonances of high order modes.

We proposed a few BHD mode-matching telescope designs and then preformed a few monte-carlo experiments to see how the imperfections would change the story. We assumed a 2 mm (1-sigma) error on the location of the components and 1% (1-sigma) fractional error on the RoC of the curved mirrors. The angle of incidence has not yet been taken into account (no astigmatism at the moment but will be included in the follow-up study.)

For the LO path things are mostly fine. We can use LO1 and LO2 as the actuators (Sec. 2.2 of the note), and when errors are taken into account more than 90% of times we can still achieve 98% mode matching. The gouy phase separation between LO1 and LO2 > 34 deg for 90% of the time, which corresponds to a condition number of the sensing matrix of ~ 3.

The situation is more tricky for the AS path. While the telescopes are usually robust against 2 or 3 mm of positional error, the 1% RoC does affect the performance quite significantly. In the note we choose two best-performing ones but still only 50% of the time they can maintain a power-overlap of > 99%. In fact, the 1% RoC error assumed should be quite optimistic... Not sure if we could achieve this in reality.

One potential way out is to ignore the MM for the first round of BHD. Here anyway we only need to test the ISC schemes. Then in the second round when we have the whole BHD board suspended, we can then use AS1 and the BHD board as the actuators. This might be able to make things more forgiving if we don't need to shrink the AS beam very fast so that it could be separated from AS4 in gouy phase.

As Rana suggested, we present the scattering plot of the AS path mode matching for various variables. The plot is for the AS path, Plan 2 (whose params we summarize at the end of this entry).

In the corner plot, we color-coded each realization according to the mode matching. We use (purple, olive, grey) for (MM>0.99, 0.98<MM<=0.99, MM<=0.98), respectively. From the plot, we can see that it is most sensitive to the RoC of AS1. The plot also shows that we can compensate for some of the MM errors if we adjust the distance between AS1-AS3 (note that AS2 is a flat mirror). The telescope is quite robust to other errors.

The compensation requirement is further shown in the second plot. To correct for the 1% RoC error of AS1, we typically need to adjust AS1-AS3 distance by ~ 1 cm (if we want to go back to MM=1; the window for >0.99 MM spans also about 1 cm). This should be doable because the nominal distance between AS1-AS3 is 115 cm.

The story for plan1 is similar and thus not shown here.

I think the conclusion is that if the AS1 RoC error is not significantly more than 1%, then with some adjustment of the AS1-AS3 distance (~ 1 cm), we could find a solution that simultaneously makes the AS path mode-matching better than 99% for the t- and s-planes.

The requirement of the LO path is less strict and the current plan using LO1-LO2 actuation should work.

We computed the required actuation range for the telescope design in elog:15357. The result is summarized in the table below. Here we assume we misalign an IFO mirror by 1 urad, and then compute how many urad do we need to move the (AS1, AS4) or (LO1, LO2) mirrors to simultaneously correct for the two gouy phases.

Actuation requirement in urad per urad misalignment

[urad/urad]

ITMX

ITMY

ETMX

ETMY

BS

PRM

PR2

PR3

SR3

SRM

AS1

1.9

2.1

-5.0

-5.5

0.5

0.5

-0.3

0.2

0.1

0.6

AS4

2.9

2.0

-8.8

-5.5

-5.9

-0.7

1.3

-0.7

-0.5

0.7

LO1

-4.0

-3.9

11.0

10.4

1.9

-0.4

-0.2

0.1

0.0

-1.1

LO2

-5.0

-3.7

15.1

10.4

8.7

0.8

1.9

1.1

0.7

-1.3

The most demanding ifo mirrors are the ETMs and the BS, for every 1 urad misalignment the telescope needs to move 10-15 urad to correct for that. However, it is unlikely for those mirrors to move more 100 nrad for a locked ifo with ASC engaged. Thus a few urad actuation should be sufficient. For the recycling mirrors, every 1 urad misalignment also requires ~ 1 urad actuation.

As a result, if we could afford 10 urad actuation range for each telescope suspension, then the gouy phase separations we have should be fine.

We looked at the oplev spectra from gps 1274418500 for 512 sec. This should be a period when the ifo was locked in the PRFPMI state according to elog:15348. We just focused on the yaw data for now. Please see the attached plots. The solid traces are for the ASD, and the dotted ones are the cumulative rms. The total rms for each mirror is also shown in the legend.

I am now confused... The ITMs looked somewhat reasonable in that at least the < 1 Hz motion was suppressed. The total rms is ~ 0.1 urad, which was what I would expect naively (~ x100 times worse than aLIGO).

There seems to be no low-freq suppression on the ETMs though... Is there no arm ASC at the moment???

We consider the astigmatism effects of the stock options. The conclusions are:

1. For the AS path, the stock should work fine for the phase-one of BHD, if we could tolerate a few percent MM loss. The window for length adjustment to achieve >99% MM for both s and t is only 1 mm for 1% RoC error (compared to ~ 1 cm in the customized case).

2. The LO path seemed tricky. As LO3 & LO4 are both significantly curved (RoC<=0.5 m), the non-zero angle of incidence makes the astigmatism quite sever. For the t-plane the nominal MM can be 0.98, yet for the s-plane, the nominal MM is only 0.72. We could move things around to achieve a MM ~ 0.85, which is probably fine for the phase-one implementation but not long term.

Details:

Attachments 1-3 are for the AS path; 4-6 are for the LO path.

1 & 4. Marginalized MM distribution for the AS/LO paths. Here we assumed 5 mm positional error and 1% fractional RoC error. Due to the astigmatism, the nominal s-plane MM is only 0.72 for the LO path.

2 & 5. Scattering plots for the AS/LO paths. We color coded the points as the following: pink: MM>0.99; olive: 0.98<MM<=0.99; grey: MM<=0.98. For the AS path, MM is mostly sensitive to the AS1 RoC and can be adjusted by changing AS1-AS3 distance. For the LO path, the LO3 RoC and LO3-LO4 distance are most critical for the MM.

3 & 6. Assuming +- 1% AS1 (LO3) fractional RoC error, how much can we compensate for it using AS1-AS3 (LO3-LO4) distance. For the AS path, there exists a ~ 1 mm window where the MM for s and t can simultaneously > 99%. For the LO path, the best we can do is to make s and t both ~ 85%.

Quote:

Summary

For the initial phase of BHD testing, we recently discussed whether the mode-matching telescopes could be built with 100% stock optics. This would allow the optical system to be assembled more quickly and cheaply at a stage when having ultra-low loss and scattering is less important. I've looked into this possibility and conclude that, yes, we do have a good stock optics option. It in fact achieves comprable performance to our optimized custom-curvature design [ELOG 15357]. I think it is certainly sufficient for the initial phase of BHD testing.

Vendor

It turns out our usual suppliers (e.g., CVI, Edmunds) do not have enough stock options to meet our requirements. This is for two reasons:

For sufficient LO1-LO2 (AS1-AS4) Gouy phase separation, we require a very particular ROC range for LO1 (AS1) of 5-6 m (2-3 m).

We also require a 2" diameter for the suspended optics, which is a larger size than most vendors stock for curved reflectors (for example, CVI has no stock 2" options).

However I found that Lambda Research Optics carries 1" and 2" super-polished mirror blanks in an impressive variety of stock curvatures. Even more, they're polished to comprable tolerances as I had specificied for the custom low-scatter optics [DCC E2000296]: irregularity < λ/10 PV, 10-5 scratch-dig, ROC tolerance ±0.5%. They can be coated in-house for 1064 nm to our specifications.

From modeling Lambda's stock curvature options, I find it still possible to achieve mode-matching of 99.9% for the AS beam and 98.6% for the LO beam, if the optics are allowed to move ±1" from their current positions. The sensitivity to the optic positions is slightly increased compared to the custom-curvature design (but by < 1.5x). I have not run the stock designs through Hang's full MC corner-plot analysis which also perturbs the ROCs [ELOG 15339]. However for the early BHD testing, the sensitivity is secondary to the goal of having a quick, cheap implementation.

Stock-Part Telescope Designs

The following tables show the best telescope designs using stock curvature options. It assumes the optics are free to move ±1" from their current positions. For comparison, the values from the custom-curvature design are also given in parentheses.

AS Path

The AS relay path is shown in Attachment 1:

AS1-AS4 Gouy phase separation: 71°

Mode-matching to OMC: 99.9%

Optic

ROC (m)

Distance from SRM AR (m)

AS1

2.00 (2.80)

0.727 (0.719)

AS2

Flat (Flat)

1.260 (1.260)

AS3

0.20 (-2.00)

1.864 (1.866)

AS4

0.75 (0.60)

2.578 (2.582)

LO Path

The LO relay path is shown in Attachment 2:

LO1-LO2 Gouy phase separation: 67°

Mode-matching to OMC: 98.6%

Optic

ROC (m)

Distance from PR2 AR (m)

LO1

5.00 (6.00)

0.423 (0.403)

LO2

1000 (1000)

2.984 (2.984)

LO3

0.50 (0.75)

4.546 (4.596)

LO4

0.15 (-0.45)

4.912 (4.888)

Ordering Information

I've created a new tab in the BHD procurement spreadsheet ("Stock MM Optics Option") listing the part numbers for the above telescope designs, as well as their fabrication tolerances. The total cost is $2.8k + the cost of the coatings (I'm awaiting a quote from Lambda for the coatings). The good news is that all the curved substrates will receive the same HR/AR coatings, so I believe they can all be done in a single coating run.

We explore bilinear SRCL to DARM noise coupling mechanisms, and show two cases that by doing BHD readout the noise performance can be improved. In the first case, the bilinear piece is due to residual DHARD motion (see also LHO:45823), and it matters mostly for the low-frequency (<100 Hz) part, and in the second piece the bilinear piece is due to residual SRCL fluctuation and it matters mostly for the a few x 100 Hz part. Details are below:

=================================================

General Model:

We can write the SRCL to DARM transfer function as (Evan Hall's thesis, eq. 2.29)

The first term in (2) is due to residual DARM motion dx_D. This term does not depends on the DC value of DARM offset <x_D> and thus does not depend on doing BHD or DC readout. On the other hand, the typical residual DARM motion is 1 fm << 1 pm of DARM offset. Since the current feedforward reduction factor is about 10 (see both Den Martynov's thesis and Evan Hall's thesis), clearly we are not limited by the residual DARM motion.

The second term is due to the change in the arm finesse, which can be affected by, e.g., the alignment fluctuation (both increasing the loss due to scattering into 01/10 modes and affecting the spot positon and hence changing the losses), and is likely to be the reason why we see the effect being modulated by DHARD.

The last term in (2) is due to the residual SRCL fluctuation and is important for the ~ a few x 100 Hz band.

=================================================

DHARD effects.

As argued above, the DHARD affects the SRCL -> DARM coupling as it changes the finesse in the arm cavity (through scattering into 01/10 modes; in finesse we cannot directly simulate the effects due to spot hitting a rougher location).

Since in the second term of eq. (2) the LF part depends on the DARM DC offset <x_D>, this effect can be improved by going from DC readout to BHD.

To simulate it in finesse, at a fixed DARM DC offset, we compute the SRCL->DARM transfer functions at different DHARD offsets, and then numerically compute the derivative \partial Z_s2d / \partial \theta_{DH}. Then multiplying this derivative with the rms value of DHARD fluctuation \theta_{DH} we then know the expected bilinear coupling piece.

The result is shown in the first attached plot. Here we have assumed a flat SRCL noise of 5e-16 m/rtHz for simplicity (see PRD 93, 112004, 2016). We do not account for the loop effects which further reduces the high frequency components for now. The residual DHARD RMS is assumed to be 1 nrad.

In the first plot, from top to bottom we show the SRCL noise projection at different DARM DC offsets of (0.1, 1, 10) pm. Since the DHARD alignment only affects the arm finesse starting at quadratic order, it thus matters what DC offset in DHARD we assume. In each pannel, the blue trace is for no DC offset in DHARD and the orange one for a 5 nrad DC offset. As a reference, the A+ sensitivity is shown in grey trace in each plot as a reference.

We can see if there is a large DC offset in DHARD (a few nrad) and we still do DC readout with a few pm of DARM offset, then the bilinear piece of SRCL can still contaminate the sensitivity in the 10-100 Hz band (bottom panel; orange trace). On the other hand, if we do BHD, then the SRCL noise should be down by ~ x100 even compared to with the top panel.

(A 5 nrad of DC offset in DHARD coupled with 1 nrad RMS would cause about 0.5% RIN in the arms. This is somewhat greater than the typically measured RIN which is more like <~ 0.2%. See the second plot).

=================================================

SRCL effect.

Similarly we can consider the SRCL->DARM coupling due to residual SRCL rms. The approach is very similar to what we did above for DHARD. I.e., we compute Z_s2d at fixed DARM offset and for different SRCL offsets, then we numerically evaluate \partial Z_s2d / \partial dphi_S. A residual SRCL rms of 0.1 nm is then used to generate the projection shown in the third figure.

Unlike the DHARD effect, the bilinear SRCL piece does not depend on the DC SRCL detuning (for the 50-500 Hz part). It does still depends on the DARM DC offset and therefore could be improved by BHD.

Since we do not include the LP of the SRCL loop in this plot, the HF noise at 1 kHz is artifical as it can be easily filtered out. However, the LP will not be very strong around 100-300 Hz for a SRCL UGF ~ 30 Hz, and thus doing BHD could still have some small improvements for this effect.

What: Anchal and I measured the XARM OLTF last Thursday.

Goal: 1. measure the 2 zeros and 2 poles in the analog whitening filter, and potentially constrain the cavity pole and an overall gain.

2. Compare the parameter distribution obtained from measurements and that estimated analytically from the Fisher matrix calculation.

3. Obtain the optimized excitation spectrum for future measurements.

How: we inject at C1:SUS-ETMX_LSC_EXC so that each digital count should be directly proportional to the force applied to the suspension. We read out the signal at C1:SUS-ETMX_LSC_OUT_DQ. We use an approximately white excitation in the 50-300 Hz band, and intentionally choose the coherence to be only slightly above 0.9 so that we can get some statistical error to be compared with the Fisher matrix's prediction. For each measurement, we use a bandwidth of 0.25 Hz and 10 averages (no overlapping between adjacent segments).

The 2 zeros and 2 poles in the analog whitening filter and an overall gain are treated as free parameters to be fitted, while the rest are taken from the model by Anchal and Paco (elog:16363). The optical response of the arm cavity seems missing in that model, and thus we additionally include a real pole (for the cavity pole) in the model we fit. Thus in total, our model has 6 free parameters, 2 zeros, 3 poles, and 1 overall gain.

The analysis codes are pushed to the 40m/sysID repo.

Fig. 1 shows one measurement. The gray trace is the data and the olive one is the maximum likelihood estimation. The uncertainty for each frequency bin is shown in the shaded region. Note that the SNR is related to the coherence as

SNR^2 = [coherence / (1-coherence)] * (# of average),

and for a complex TF written as G = A * exp[1j*Phi], one can show the uncertainty is given by

\Delta A / A = 1/SNR, \Delta \Phi = 1/SNR [rad].

Fig. 2. The gray contours show the 1- and 2-sigma levels of the model parameters using the Fisher matrix calculation. We repeated the measurement shown in Fig. 1 three times, and the best-fit parameters for each measurement are indicated in the red-crosses. Although we only did a small number of experiments, the amount of scattering is consistent with the Fisher matrix's prediction, giving us some confidence in our analytical calculation.

One thing to note though is that in order to fit the measured data, we would need an additional pole at around 1,500 Hz. This seems a bit low for the cavity pole frequency. For aLIGO w/ 4km arms, the single-arm pole is about 40-50 Hz. The arm is 100 times shorter here and I would naively expect the cavity pole to be at 3k-4k Hz if the test masses are similar.

Fig. 3. We then follow the algorithm outlined in Pintelon & Schoukens, sec. 5.4.2.2, to calculate how we should change the excitation spectrum. Note that here we are fixing the rms of the force applied to the suspension constant.

Fig. 4 then shows how the expected error changes as we optimize the excitation. It seems in this case a white-ish excitation is already decent (as the TF itself is quite flat in the range of interest), and we only get some mild improvement as we iterate the excitation spectra (note we use the color gray, olive, and purple for the results after the 0th, 1st, and 2nd iteration; same color-coding as in Fig. 3).

Yesterday afternoon Paco and I measured the PRM L2P transfer function. We drove C1:SUS-PRM_LSC_EXC with a white noise in the 0-10 Hz band (effectively a white, longitudinal force applied to the suspension) and read out the pitch response in C1:SUS-PRM_OL_PIT_OUT. The local damping was left on during the measurement. Each FFT segment in our measurement is 32 sec and we used 8 non-overlapping segments for each measurement. The empirically determined results are also compared with the Fisher matrix estimation (similar to elog:16373).

Results:

Fig. 1 shows one example of the measured L2P transfer function. The gray traces are measurement data and shaded region the corresponding uncertainty. The olive trace is the best fit model.

Note that for a single-stage suspension, the ideal L2P TF should have two zeros at DC and two pairs of complex poles for the length and pitch resonances, respectively. We found the two resonances at around 1 Hz from the fitting as expected. However, the zeros were not at DC as the ideal, theoretical model suggested. Instead, we found a pair of right-half plane zeros in order to explain the measurement results. If we cast such a pair of right-half plane zeros into (f, Q) pair, it means a negative value of Q. This means the system does not have the minimum phase delay and suggests some dirty cross-coupling exists, which might not be surprising.

Fig. 2 compares the distribution of the fitting results for 4 different measurements (4 red crosses) and the analytical error estimation obtained using the Fisher matrix (the gray contours; the inner one is the 1-sigma region and the outer one the 3-sigma region). The Fisher matrix appears to underestimate the scattering from this experiment, yet it does capture the correlation between different parameters (the frequencies and quality factors of the two resonances).

One caveat though is that the fitting routine is not especially robust. We used the vectfit routine w/ human intervening to get some initial guesses of the model. We then used a standard scipy least-sq routine to find the maximal likelihood estimator of the restricted model (with fixed number of zeros and poles; here 2 complex zeros and 4 complex poles). The initial guess for the scipy routine was obtained from the vectfit model.

Fig. 3 shows how we may shape our excitation PSD to maximize the Fisher information while keeping the RMS force applied to the PRM suspension fixed. In this case the result is very intuitive. We simply concentrate our drive around the resonance at ~ 1 Hz, focusing on locations where we initially have good SNR. So at least code is not suggesting something crazy...

Fig. 4 then shows how the new uncertainty (3-sigma contours) should change as we optimize our excitation. Basically one iteration (from gray to olive) is sufficient here.

We will find a time very recently to repeat the measurement with the optimized injection spectrum.

We tried to compare the parameter estimation uncertainties of white vs. optimal excitation. We drove C1:SUS-PRM_LSC_EXC with "Normal" excitation and digital gain of 700.

For the white noise exciation, we simply put a butter("LowPass",4,10) filter to select out the <10 Hz band.

For the optimal exciation, we use butter("BandPass",6,0.3,1.6) gain(3) notch(1,20,8) to approximate the spectral shape reported in elog:16384. We tried to use awg.ArbitraryLoop yet this function seems to have some bugs and didn't run correctly; an issue has been submitted to the gitlab repo with more details. We also noticed that in elog:16384, the pitch motion should be read out from C1:SUS-PRM_OL_PIT_IN1 instead of the OUT channel, as there are some extra filters between IN1 and OUT. Consequently, the exact optimal exciation should be revisited, yet we think the main result should not be altered significantly.

While a more detail analysis will be done later offline, we post in the attached plot a comparison between the white (blue) vs optimal (red) excitation. Note in this case, we kept the total force applied to the PRM the same (as the RMS level matches).

Under this simple case, the optimal excitation appears reasonable in two folds.

First, the optimization tries to concentrate the power around the resonance. We would naturally expect that near the resonance, we would get more Fisher information, as the phase changes the fastest there (i.e., large derivatives in the TF).

Second, while we move the power in the >2 Hz band to the 0.3-2 Hz band, from the coherence plot we see that we don't lose any information in the > 2 Hz region. Indeed, even with the original white excitation, the coherence is low and the > 2 Hz region would not be informative. Therefore, it seems reasonable to give up this band so that we can gain more information from locations where we have meaningful coherence.

We report here the analysis results for the measurements done in elog:16388.

Figs. 1 & 2 are respectively measurements of the white noise excitation and the optimized excitation. The shaded region corresponds to the 1-sigma uncertainty at each frequency bin. By eyes, one can already see that the constraints on the phase in the 0.6-1 Hz band are much tighter in the optimized case than in the white noise case.

We found the transfer function was best described by two real poles + one pair of complex poles (i.e., resonance) + one pair of complex zeros in the right-half plane (non-minimum phase delay). The measurement in fact suggested a right-hand pole somewhere between 0.05-0.1 Hz which cannot be right. For now, I just manually flipped the sign of this lowest frequency pole to the left-hand side. However, this introduced some systematic deviation in the phase in the 0.3-0.5 Hz band where our coherence was still good. Therefore, a caveat is that our model with 7 free parameters (4 poles + 2 zeros + 1 gain as one would expect for an ideal signal-stage L2P TF) might not sufficiently capture the entire physics.

In Fig. 3 we showed the comparison of the two sets of measurements together with the predictions based on the Fisher matrix. Here the color gray is for the white-noise excitation and olive is for the optimized excitation. The solid and dotted contours are respectively the 1-sigma and 3-sigma regions from the Fisher calculation, and crosses are maximum likelihood estimations of each measurement (though the scipy optimizer might not find the true maximum).

Note that the mean values don't match in the two sets of measurements, suggesting potential bias or other systematics exists in the current measurement. Moreover, there could be multiple local maxima in the likelihood in this high-D parameter space (not surprising). For example, one could reduce the resonant Q but enhance the overall gain to keep the shoulder of a resonance having the same amplitude. However, this correlation is not explicit in the Fisher matrix (first-order derivatives of the TF, i.e., local gradients) as it does not show up in the error ellipse.

In Fig. 4 we show the further optimized excitation for the next round of measurements. Here the cyan and olive traces are obtained assuming different values of the "true" physical parameter, yet the overall shapes of the two are quite similar, and are close to the optimized excitation spectrum we already used in elog:16388.

We did a few quick XARM oltf measurements. We excited C1:LSC-ETMX_EXC with a broadband white noise upto 4 kHz. The timestamps for the measurements are: 1318199043 (start) - 1318199427 (end).

We will process the measurement to compute the cavity pole and analog filter poles & zeros later.

One goal of our sysID study is to improve the aLIGO L2A feedforward. Our algorithm currently improves only the statistical uncertainty and assumes the systematic errors are negligible. However, I am currently baffled by how to fit a (nearly) realistic suspension model...

My test study uses the damped aLIGO QUAD suspension model. From the Matlab model I extract the L2 drive in [N] to L3 pitch in [rad] transfer function (given by a SS model with the A matrix having a shape of 103x103). I then tried to use VectFIT to fit the noiseless TF. After removing nearby z-p pairs (defined by less than 0.2 times the lowest pole frequency) and high-frequency zeros, I got a model with 6 complex pole pairs and 4 complex zero pairs (21 free parameters in total). I also tried to fit the TF (again, noiseless) with an MCMC algorithm assuming the underlying model has the same number of parameters as the VectFIT results.

Please see the first attached plots for a comparison between the fitted models and the true one. In the second plot, we show the fractional residual

| TF_true - TF_fit | / | TF_true |,

and the inverse of this number gives the saturating SNR at each frequency. I.e., when the statistical SNR is more than the saturating value, we are then limited by systematic errors in the fitting. And so far, disappointingly I can only get an SNR of 10ish for the main resonances...

I wonder if people know better ways to reduce this fitting systematic... Help is greatly appreciated!

We have been discussing how does the parameter estimation depends on the length per FFT segment. In other words, after we collected a series of data, would it be better for us to divide it into many segments so that we have many averages, or should we use long FFT segments so that we have more frequency bins?

My conclusions are that:

1). We need to make sure that the segment length is long enough with T_seg > min[ Q_i / f_i ], where f_i is the resonant frequency of the i'th resonant peak and the Q_i its quality factor.

2). Once 1) is satisfied, the result depends weakly on the FFT length. There might be a weak hint preferring a longer segment length (i.e., want more freq bins than more averages) though.

To reach the conclusion, I performed the following numerical experiment.

I considered a simple pendulum with resonant frequency f_1 = 0.993 Hz and Q_1 = 6.23. The value of f_1 is chosen such that it is not too special to fall into a single freq bin. Additionally, I set an overall gain of k=20. I generated T_tot = 512 s of data in the time domain and then did the standard frequency domain TF estimation. I.e., I computed the CSD between excitation and response (with noise) over the PSD of the excitation. The spectra of excitation and noise in the readout channel are shown in the first plot.

In the second plot, I showed the 1-sigma errors from the Fisher matrix calculation of the three parameters in this problem, as well as the determinant of the error matrix \Sigma = inv(Fisher matrix). All quantities are plotted as functions of the duration per FFT segment T_seg. The red dotted line is [Q_1/f_1], i.e., the time required to resolve the resonant peak. As one would expect, if T_seg <~ (Q_1/f_1), we cannot resolve the dynamics of the system and therefore we get nonsense PE results. However, once T_seg > (Q_1/f_1), the PE results seem to be just fluctuating (as f_1 does not fall exactly into a single bin). Maybe there is a small hint that longer T_seg is better. Potentially, this might be due to that we lose less information due to windowing? To be investigated further...

I also showed the Fisher estimation vs. MCMC results in the last two plots. Here each dot is an MCMC posterior. The red crosses are the true values, and the purple contours are the results of the Fisher calculations (3-sigma contours). The MCMC results showed similar trends as the Fisher predictions and the results for T_seg = (32, 64, 128) s all have similar amounts of scattering << the scattering of the T_seg=8 s results. Though somehow it showed a biased result. In the third plot, I manually corrected the mean so that we could just compare the scattering. The fourth plot showed the original posterior distribution.

1. In the error propogation equation, it should be \Delta \Theta = -H^{-1} M \Delta \Lambda, instead of the fractional error.

2. For the astro parameters, in general you would need t_c for the time of coalescence and \phi_c for the phase. See, e.g., https://ui.adsabs.harvard.edu/abs/1994PhRvD..49.2658C/abstract.

3. Fig. 1 looks very nice to me, yet I don't understand Fig. 3... Why would phase or amplitude uncertainties at 30 Hz affect the tidal deformability? The tide should be visible only > 500 Hz.

4. For BBH, we don't measure individual spin well but only their mass-weighted sum, \chi_eff = (m_1*a_1 + m_2*a_2)/(m_1 + m_2). If you treat S1z and S2z as free parameters, your matrix is likely degenerate. Might want to double-check. Also, for a BBH, you don't need to extend the signal much higher than \omega ~ 0.4/M_tot ~ 10^4 Hz * (Ms/M_tot). So if the total mass is ~ 100 Ms, then the highest frequency should be ~ 100 Hz. Above this number there is no signal.

I modified astroFisherLib.py to include these parameters. Please note that their meaning is that we don't know when the signal happens and at which phase it merges.

It does not mean the time & phase from a reference frequency to the merger. This part is not free to vary because it is fixed by the intrinsic parameters.

It might be good to have a quick scan through the Cutler & Flanagan 94 paper to better understand their physical meanings.

To use a razorblade to measure beam waist at multiple points along the optical axis, so as to later extrapolate the modal profile of the entire beam. This information will then be used to effectively couple AUX laser light to fibers for use in the frequency offset locking apparatus.

Data Acquisition

1) Step the micrometer-controlled razorblade across the beam at a given value of Z, along optical axis, in the plane orthogonal to it (arbitrarily called X).

2) At each value of X, record the corresponding output of a photodiode, (Thorlabs PD A55) here given in mV.

3) Repeat process at multiple points along Z

Analysis

Data from each iteration in the X were fitted to the error function shown below.

V(x) = A*(erf((x-m)/s)+c)

In the Y, they were fitted to:

V(x) = -A*(erf((x-m)/s)+c)

'A' corresponds to an amplitude, 'm' to a mean, 's' to a σ, and 'c' to an offset.

(Only because in Y measurements, the blade progressed toward eclipsing the beam, as opposed to in the X where it progressively revealed the beam.

These fits can be solved for x = (erf^{-1}((V/A)-c)*s)+m1 which can be calculated at the points (V_{max}/e^{2}) and (V_{max}*(1-1/e^{2})). The difference between these points will yield beam waist, w(z).