Graduate Studies, College of (Okanagan)
DSpace
UBCO
Yeremi, Miayan
Master of Science - MSc
University of British Columbia
We have developed an experimentally designed statistical framework, which allows for observation and simulation position-position-velocity cubes to be compared in a systematic manner. The comparison between cubes can be performed using a multitude of metric measures, in particular we have used the following metrics: principal component analysis (PCA), spectral correlation function (SCF), Cramer, and maximum mean discrepancy (MMD). We ran a suite of carefully chosen simulations, and performed a simulation-to-simulation comparison to test the framework in a controlled environment. We have tested the effects of four physical parameters, Mach number, magnetic field, driving scale, and temperature. The developed framework extracted key physics from the optimally designed experiment. In addition, merits of the different measures were highlighted in a standardized manner. After all four metrics were analyzed, we have concluded that PCA, Cramer and MMD respond in a linear manner to changes in factor settings, where as SCF displayed signs of curvature. Furthermore, all four statistics picked out a consistent physical trend, predicting that Mach num- ber, temperature, and Mach-number-temperature interaction are significant dissimilarity drivers. Finally, we propose that this framework will be used in future work for the calibration of simulations to observational data.
Statistically Calibrating ComparisonMetrics of Star Formation SimulationsbyMiayan YeremiB.Sc., The University of British Columbia, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate Studies(IGS)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)July 2013c? Miayan Yeremi, 2013AbstractWe have developed an experimentally designed statistical framework,which allows for observation and simulation position-position-velocity cubesto be compared in a systematic manner. The comparison between cubescan be performed using a multitude of metric measures, in particular wehave used the following metrics: principal component analysis (PCA), spec-tral correlation function (SCF), Cramer, and maximum mean discrepancy(MMD). We ran a suite of carefully chosen simulations, and performed asimulation-to-simulation comparison to test the framework in a controlledenvironment. We have tested the effects of four physical parameters, Machnumber, magnetic field, driving scale, and temperature. The developedframework extracted key physics from the optimally designed experiment.In addition, merits of the different measures were highlighted in a standard-ized manner. After all four metrics were analyzed, we have concluded thatPCA, Cramer and MMD respond in a linear manner to changes in factorsettings, where as SCF displayed signs of curvature. Furthermore, all fourstatistics picked out a consistent physical trend, predicting that Mach num-ber, temperature, and Mach-number-temperature interaction are significantdissimilarity drivers. Finally, we propose that this framework will be usedin future work for the calibration of simulations to observational data.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . ixChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 1Chapter 2: Theoretical Background . . . . . . . . . . . . . . . . 92.1 The Physics of Star Formation . . . . . . . . . . . . . . . . . 92.2 Enzo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.3 Post-Processing . . . . . . . . . . . . . . . . . . . . . . . . . . 162.4 Comparison Statistic . . . . . . . . . . . . . . . . . . . . . . . 182.4.1 Principal Component Analysis . . . . . . . . . . . . . 192.4.2 The Spectral Correlation Function . . . . . . . . . . . 202.4.3 Cramer statistic . . . . . . . . . . . . . . . . . . . . . 202.4.4 Maximum Mean Discrepancy . . . . . . . . . . . . . . 212.5 Experimental Design . . . . . . . . . . . . . . . . . . . . . . . 212.6 Response Analysis . . . . . . . . . . . . . . . . . . . . . . . . 22Chapter 3: Methodology . . . . . . . . . . . . . . . . . . . . . . 233.1 Physical Parameters/Factors ran in the simulation . . . . . . 233.1.1 Mach Number . . . . . . . . . . . . . . . . . . . . . . 233.1.2 Magnetic Field . . . . . . . . . . . . . . . . . . . . . . 243.1.3 Driving scale . . . . . . . . . . . . . . . . . . . . . . . 243.1.4 Temperature . . . . . . . . . . . . . . . . . . . . . . . 253.2 Physical Fixed Parameters ran in the simulation . . . . . . . 25iiiTABLE OF CONTENTS3.3 Data Collection and Analysis Summary . . . . . . . . . . . . 263.4 Experimental Designs . . . . . . . . . . . . . . . . . . . . . . 29Chapter 4: Results . . . . . . . . . . . . . . . . . . . . . . . . . . 304.1 Baseline Design . . . . . . . . . . . . . . . . . . . . . . . . . 304.2 Orthogonal Design . . . . . . . . . . . . . . . . . . . . . . . . 374.2.1 Fractional Factorial Orthogonal Design . . . . . . . . 384.3 Unreplicated Design Analysis . . . . . . . . . . . . . . . . . . 394.3.1 Normal Plots/Half-Normal Plots . . . . . . . . . . . . 394.3.2 Centre Point Variance Estimation . . . . . . . . . . . 404.3.3 Lenth?s Method . . . . . . . . . . . . . . . . . . . . . 43Chapter 5: Conclusion/Discussion . . . . . . . . . . . . . . . . . 495.1 Result Summary . . . . . . . . . . . . . . . . . . . . . . . . . 495.2 Result Interpretation . . . . . . . . . . . . . . . . . . . . . . . 505.3 Future investigation of Curvature . . . . . . . . . . . . . . . . 515.4 Future Research Directions . . . . . . . . . . . . . . . . . . . 51Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54ivList of TablesTable 3.1 The parameters listed in this table were held fixedacross simulation runs. . . . . . . . . . . . . . . . . . . 26Table 3.2 Cramer statistic comparing a fiducial run?s specifictime step to all the rest of the 15 simulations. Zero iscoded for the low setting of each parameter, and oneis conversely coded for the high setting. . . . . . . . . 27Table 4.1 Baseline factor settings . . . . . . . . . . . . . . . . . . 30Table 4.2 One factor at a time baseline design results for PCAstatistic. The first column gives the effect estimatesand the last column gives the p-value, or measure ofsignificance. The smaller the p-value the more signifi-cant the corresponding effect estimate. . . . . . . . . . 32Table 4.3 PCA full factorial saturated model linear regressionresults. . . . . . . . . . . . . . . . . . . . . . . . . . . . 34Table 4.4 Fractional factorial for PCA statistic under baselineparameterization regression results. . . . . . . . . . . . 35Table 4.5 Lack of Fit p-values for the four different statistics, ordistance metrics. Smaller p-values mean that there ismore evidence for curvature. Note, p-values are prob-abilities. Therefore, they range between 0 and 1. Theonly metric that might exhibit signs of curvature is SCF. 36Table 4.6 Significant effects in Fractional Factorial Design. . . . 40Table 4.7 Significant effects in Full Factorial Design. . . . . . . . 43Table 4.8 Active Factors in Fractional Factorial Design usingLenth?s Method. . . . . . . . . . . . . . . . . . . . . . 48Table 4.9 Active Factors in Full Factorial Design using Lenth?sMethod. . . . . . . . . . . . . . . . . . . . . . . . . . . 48vList of FiguresFigure 1.1 Star Formation Process. This image was found andproduced in ref. [16]. . . . . . . . . . . . . . . . . . . 2Figure 1.2 Hubble telescope image of the Horsehead GMC. Thisimage was found and produced in ref. [26]. . . . . . . 3Figure 1.3 Density simulation slices through the 3D position-position-position (PPP) space of the simulations, be-fore they are converted to position-position-velocity(PPV) data. The four panels are taken from twodifferent runs with different high and low initial mag-netic field setting. Each run has two panels, whichcorrespond to different orientations of slicing the 3DPPP space. This slices show that the magnetic fieldparameter has a spatial directionality dependence. . . 6Figure 2.1 The top three plots are of the Multi-Grid (MG) methodfor three different iterations. The bottom three plotsare of the Single-Grid (SG) method for three differentiterations. The last plot for each respective methodhas roughly the same measure of error, but the MGmethod has taken fewer iterations. . . . . . . . . . . 15viLIST OF FIGURESFigure 3.1 Examples of PPV slices extracted from four differentsimulation runs. Panel (a) is the fiducial run whichall comparisons are made back to. Fiducial run hasMach number of 6, magnetic field of 1 ?G, defaultwider driving scale, and a temperature of 10 K. Panel(b) has Mach number set to 18 and all the rest of theparameters as seen in the fiducial run. Panel (c) hasthe smaller driving scale, which corresponds to thehigh setting of driving scale parameter. Lastly, panel(d) shows a slice for high magnetic field of 100 ?Gand all the rest of the parameters are set at theirrespective fiducial settings. This image was generatedin ref. [28]. . . . . . . . . . . . . . . . . . . . . . . . . 28Figure 4.1 Several factors that interact across the different dis-tance metrics. If the two lines are not parallel then aninteraction between the labeled parameters is present. 33Figure 4.2 Effect estimates half normal plot for fractional facto-rial model. The dots that are sitting high above therest or effectively outliers are considered as significanteffects. I have labeled the two most significant effects. 41Figure 4.3 Effect estimates half normal plot for the full facto-rial model. The dots that are sitting high above therest or effectively outliers are considered as significanteffects. I have labeled the two most significant effects. 42Figure 4.4 PCA and SCF Lenth results for fractional factorialdesign. Note, that any bar protruding out of the twohorizontal lines is considered significant. The signifi-cant main effect for PCA is Mach and the significantmain effects for SCF are Mach and Temp. PCA haslost the temperature significance moving between thefull and fractional factorial, see Figure 4.6 to high-light the difference. . . . . . . . . . . . . . . . . . . . 44viiLIST OF FIGURESFigure 4.5 Cramer and MMD Lenth results for fractional fac-torial design. Significance is judged as described inprevious caption. Note, the significant main effectfor MMD is Mach and the significant main effectsfor Cramer are Mach and Temp. MMD has lost thetemperature significance moving between the full andfractional factorial, see Figure 4.7 to highlight thedifference. . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 4.6 PCA and SCF Lenth results for full factorial design.Note, that any bar protruding out of the two hori-zontal lines is considered significant. The significantmain effects for both PCA and SCF are Mach andTemp. I omitted the higher order interactions tomake the plot more readable; they were not signif-icant. . . . . . . . . . . . . . . . . . . . . . . . . . . . 46Figure 4.7 Cramer and MMD Lenth results for full factorial de-sign. Significance is judged as described in previouscaption. Note, the main effects that are significantfor Cramer and MMD are the same as for PCA andSCF. I omitted the higher order interactions to makethe plot more readable; they were not significant. . . 47Figure 5.1 Example of higher order curvature. The test I ranfor curvature, tests if the centre point falls on a linedrawn through the edges of the parameter space. Thisis a good test for second order curvature, but if thesurface actually exhibits third order curvature thetest might fail at detecting curvature. This can bebetter understood from the above plot. . . . . . . . . 52viiiAcknowledgementsI would like to recognize the people who have made my research possible.First, I would like to start by thanking my two awesome supervisors, JasonLoeppky and Erik Rosolowsky. They have given me the opportunity towork on a problem that is vast and extremely challenging. They have putup with my consistent pestering, and have guided me through the processof research. I have gained a huge amount of knowledge from them on ainterdisciplinary subject. Next, I would like to acknowledge my peers inthe AstroStat group, Tyler Berg and Mallory J Flynn. Tyler has helpedmeet deadlines on several occasions. Mallory has been an indispensableresearch colleague and added a unique perspective. Lastly, I would liketo mention our external collaboration with Stella Offner, who has helpedsetup the numerical simulation component of our project, and has guidedour intuition with Enzo.ixChapter 1IntroductionFun only begins when one knows the rules. In the case of ouruniverse we are currently still in progress of reading the instruc-tion manual.-Richard FeynmanStar formation is a complex physical process. Figure 1.1 is a visual in-troduction to the different terminology and phases encountered in a star?sbirth process. The birth of a star and its evolution is a cyclic process - theprocess both begins and ends in a mix of gas and dust. The gas spans thespace between stars, known as the interstellar medium (ISM). The ISM isan enormous and complex environment which we use as our physical ex-perimental laboratory to study and test astrophysical theories. The keyphysical features of the ISM are turbulent gas motions, shocks, and giantgalactic magnetic fields, threaded through interstellar space. The magne-tohydrodynamic equations compactly encompass the complex behaviour ofthese features, and hold the key to further understanding of essential starformation mechanisms. In addition, the presence of dust, molecules, atoms,ions, and free electrons means we must add the consideration of radiativetransfer to the mix.Molecular clouds are embedded within the ISM. There is a continuumof different types of molecular clouds in the ISM, but a broad classificationof these clouds has been established. Lower mass clouds are called diffusemolecular clouds; these clouds are also referred to as translucent molecularclouds. Giant molecular clouds (GMC) are significantly larger and moremassive than diffuse molecular clouds. The Horsehead is a small part ofthe Orion GMC, see Figure 1.2. The GMCs exhibit fragmented structure,where dense clumps and even denser cores are found.Located outside large molecular clouds are nearly spherical clouds, knownas Bok globules. Bok Globules are cores whose diffuse material have beenablated away. Both molecular clouds and Bok globules are considered fertilegrounds for the birth of stars.1Chapter 1. IntroductionFigure 1.1: Star Formation Process. This image was found and produced inref. [16].2Chapter 1. IntroductionFigure 1.2: Hubble telescope image of the Horsehead GMC. This image wasfound and produced in ref. [26].Stars are born with a collapse of a dense gas core. I consider four im-portant factors which drive or inhibit the collapse of dense cores - gas tem-perature, self gravity, magnetic fields embedded in the weakly ionized gas,and supersonic turbulence.Theoretically, star formation is simplified by constructing an initial modelof the collapse process, which considers only gas temperature and self grav-ity. In the absence of turbulence and magnetic field, the stability of a molec-ular cloud is determined by the balance between the internal pressure of thecloud and the gravitational force exerted by the cloud on itself (i.e., self-gravity). The Jeans length, a measure of cloud stability, is(15kT4pi?GmH)0.5,where T , ?, k, G, and mH are the gas temperature, gas density, Boltzmann?sconstant, gravitational constant, and hydrogen?s molecular mass, respec-tively. Stable molecular clouds are classified by a Jeans length that is largerthan the spatial scale of the cloud, and conversely for unstable clouds. The3Chapter 1. Introductionformer Jeans length proportionality is recast to arrive at a proportionalityfor the Jeans mass, MJ ? T32 ??12 , where MJ , T , and ?, are Jeans mass, tem-perature, and density of the system, respectively. Similar to Jeans length,Jeans mass provides a stability criterion for molecular clouds. If a cloudhas a larger mass than its Jeans mass, then it will collapse; otherwise, thecloud is considered stable. In the event of a collapse, the cloud is consideredoptically thin in the initial phase of the collapse. While a cloud is opticallythin, it maintains constant temperature as excess potential energy is radi-ated away from it. During the collapse, the density of the cloud increases.Constant temperature and increasing density leads to shrinking Jeans mass,which reinforces the collapse. In theory, this suggests that a collapsing cloudcan result in the formation of one massive star. However, observations showthat this does not commonly happen. In particular, we observe that themajority of stars formed have moderate masses on the order of our sun, andtypically form in groups, resulting in binary systems and large clusters. Thissuggests that sometime during cloud collapse, fragmentation takes place.The process of cloud fragmentation is researched extensively, and thecurrent theory is not complete. The above analysis supports the conceptof fragmentation in a simplistic manner. Jeans criterion, under the as-sumption of isothermality, suggests that inhomogeneities in the cloud areaccentuated and made locally unstable. Hence, theory suggests initial in-homogeneities locally collapse and fragment, and the fragments continuecollapsing and fragmenting. This hierarchical collapse does not continue in-definitely. There comes a stage where the isothermal collapse assumption isinsufficient. This occurs when the potential energy is not transported awayfrom the core efficiently, causing the collapse to transition from isothermalcollapse to an adiabatic one. The adiabatic collapse is considered when thecloud transition from optically thin to optically thick. In an adiabatic col-lapse the temperature of the cloud increases with an increase in density inthe following way T ? ???1, where ? is the adiabatic index, or the ratio ofspecific heats ( i.e., cPcV ). During the adiabatic collapse, Jeans length andJeans mass increase and the collapse is stabilized. The argument presentedis highly simplified, but it highlights some of the key features of the collapse.If angular momentum is also considered, the rotation of the cloud duringthe collapse results in a disk-like structure surrounding a dense central core,known as the protostar.The previous discussion describes a regime in which the star is feeding offof the potential energy initially stored in the molecular cloud. During thisperiod, the star is known as pre-main-sequence. Next, it will evolve into azero-age-main-sequence star, which is the transition from pre-main-sequence4Chapter 1. Introductionto the main-sequence stage. Before the main-sequence stage, stars are fuelledsolely by potential energy from the cloud, and the classification into themain-sequence is characterized by the beginning of hydrogen fusion, wherethe star will convert hydrogen into helium as a source of energy. Finally,after the star has exhausted its hydrogen supply for fusion, it moves intothe post-main-sequence, where it completes its life cycle, with a variety ofoutcomes, depending on the cores initial mass, and eventually the materialreturns to the ISM.As previously discussed, our analysis has ignored important physicalfactors that greatly complicate this basic picture. We must return to theconsideration of two of the most important factors in this process, whichare turbulence and magnetic fields. The magnetic field is created by themoving ions and electrons within the molecular cloud, and affects the tra-jectory of these charged particles. The turbulence in the cloud and thegravitational collapse act to shear and compress the magnetic field, causinga magnetic pressure proportional to ?B2 which slows down the collapsetimescale. However, the magnetic field only acts on charged particles, thusit can not directly slow the neutral particles in the collapse, which results indifferent collapse time scales between the two types of particles. The neutraland charged particles interact through the collision process, resulting in anincrease in the average timescale of the collapse, which we call ambipolardiffusion.Turbulence is a chaotic process which converts the fluid?s kinetic energyto thermal energy, by essentially ?stirring? the material in the cloud. Forexample, molecular cloud collisions and nearby supernova explosions resultin large scale turbulent motion in a cloud. This large scale kinetic energyis hierarchically transferred down to a micro-scale by division of circulationeddies, as referred to in section 3.1.3. Turbulence is a complicated process,especially in the extreme conditions encountered in star formation. It isthought to both inhibit and accelerate the star formation process. In someregions of the cloud it compresses the fluid, providing dense cores ripe forcollapse, whereas in different regions it acts as an effective pressure inhibitingthe collapse. To examine a detailed description of the gravitational collapseprocess, we must numerically solve the magnetohydrodynamical equations.The plots shown in Figure 1.3 illustrate density maps extracted from asubset of the numerical simulations ran for this dissertation. Figures 1.3cand 1.3d are density maps for high magnetic field of 10?4 G. In contrast,Figures 1.3a and 1.3b are density maps for low magnetic field of 10?6 G.The two maps for each magnetic field setting show a different directionalprojection of the density field, where the z-direction is parallel to the initial5Chapter 1. Introduction(a) x projection of density for a low mag-netic field(b) z projection of density for a low mag-netic field(c) x projection of density for a high mag-netic field(d) z projection of density for a high mag-netic fieldFigure 1.3: Density simulation slices through the 3D position-position-position (PPP) space of the simulations, before they are converted toposition-position-velocity (PPV) data. The four panels are taken from twodifferent runs with different high and low initial magnetic field setting. Eachrun has two panels, which correspond to different orientations of slicing the3D PPP space. This slices show that the magnetic field parameter has aspatial directionality dependence.6Chapter 1. Introductionmagnetic field and the x-direction is perpendicular to it. Based on a visualinspection, it is apparent that the two sets of simulations exhibit differentcloud structure.The idea of comparing astrophysical hydrodynamical simulations to ob-servations has been investigated by many. For example, ref. [19] computedthe spectral correlation maps for both real 13CO J = 1? 0 maps of molec-ular clouds and for different theoretical models; they used the observationdata to test the performance of these theoretical models. Another exampleis presented in ref. [3], where the authors ran a large hydrodynamical simu-lation of star cluster formation and found that hydrodynamical/sink particlesimulations can reproduce many of the observed stellar properties.One challenge of comparing observations to simulations stems from thenature of their representation. Observations are recorded in a PPV (i.e.,position-position-velocity) space, whereas simulation data is represented inPPP (i.e., position-position-poistion) space. Radiative transfer post pro-cessing is applied to simulation data to emulate observational data, as dis-cussed in Section 2.3. In most cases, comparisons of PPV cubes leverageunderlying statistics to construct a comparison metric. For example, a com-mon route is to find a population of objects within the simulation (e.g. densemolecular cores) and compare the properties of the population to those ofthe known population. This method is effective when a well-defined knownpopulation is available. In contrast, theoretical understanding of turbulencehas been used to construct similarity statistics. This method can be usedwith or without known population, see ref. [28] for details of the similaritystatistics, also known as distance metrics.Researchers use different distance metrics to establish similarity betweenobservations and simulations. Usually, the simulation parameters are variedin a one-factor-at-a-time (OFAT) design. The following are several exam-ples in which OFAT designs were used in the literature, ref. [24], ref. [22],ref. [1], ref. [6]. An OFAT design is not suitable for examining interactionsbetween parameters. In Chapter 4, I will demonstrate that interactions aresignificant and should not be ignored. I propose a systematic way to per-form comparisons between simulations and observations. Instead of runningOFAT designs, it is shown that an orthogonal factorial design allows forthe extraction of relevant significant parameters and their interactions. Inparticular, I investigate the effects of Mach number, magnetic field, drivingscale, and temperature on a suite of simulations. The parameters I investi-gated are a subset of physical parameters that I believed would be interestingto investigate; in Chapter 3 I describe the physical interpretation of each pa-rameter. Orthogonal factorial designs are powerful statistical methods that7Chapter 1. Introductionensure that a parameter space is searched in a systematic manner, ref. [27]describes the optimization and analysis of experimental designs and shouldbe referred to for an in depth discussion.Chapters 2 and 3 will explain the theoretical background and practicalsetup of 16 simulation runs, specifically chosen based on two-level orthogo-nal design. Next, four comparison metrics outlined in Chapter 2 are used tocompare the 16 runs to a fiducial simulation run. The comparison metricsgenerate four responses, which are ran through a statistical analysis pre-sented in Chapter 4. The analysis of Chapter 4 is used to evaluate thequality of the metrics. In conclusion, this work serves as the basis for futureefforts, which I will outline in Chapter 5.8Chapter 2Theoretical BackgroundThe main focus of my research was to design a framework that comparesnumerical simulations to observational data. The framework is composed ofseveral stages, as follows: simulation of the star formation process, post pro-cessing of simulation data into PPV cubes; finally, the comparison of eithersimulation to simulation, or simulation to observation PPV cubes, using avariety of similarity metrics. The comparison framework is followed by astatistical analysis of the predetermined experimental design, allowing forthe validation of the different metrics used, and the extraction of importantphysical factors. All stages involved have substantial theoretical backgroundthat I will introduce in this chapter.2.1 The Physics of Star FormationThe process of star formation is complicated by many factors. Thedense molecular clouds, which are home to stellar nurseries, are set in ex-treme physical conditions. These harsh conditions lead to physical equationsmodelling complex behaviours such as, self gravitation, chaotic turbulence,gaseous supersonic magnetized plasma and more, see [13] for more details.In this thesis, I used the magneto-hydrodynamic (MHD) equations to modelthe process of star formation. The MHD equations consist of the followingset of equations:The continuity equation,???t+? ? ?u = 0,where ? and u are the mass density and fluid?s bulk velocity. The continuityequation ensures that the mass is conserved.The momentum equation,?DuDt= ??2u+1?0(B ??)B? (12?0?B2 +?P )? ???,where t, ?, ?0, B, P , and ? are time, coefficient of viscosity, vacuum perme-ability, magnetic vector field, pressure scalar field, and gravitational scalar92.1. The Physics of Star Formationpotential field, respectively. The momentum equation ensures that momen-tum is conserved. The DDt =??t + u ? ? is the material derivative of thebulk velocity. The term u ??u is non-linear and is responsible for significantcomplications in fluid dynamics (i.e., chaotic turbulence). The negativesign in front of the pressure gradient can be explained by noting that fluidflows from high pressure to low pressure regions ( a similar argument can bemade for gravitation potential term ). The B ??B term is a stress impartedon the fluid due to the curvature in the magnetic field, also known as mag-netic tension. The ?B2 is the gradient of isotropic magnetic pressure actingon the fluid. Both magnetic terms are a consequence of the Lorentz force,assuming zero resistivity and no Hall effect (i.e., electric fields in the fluid?sframe of reference are negligible). Lastly, ??? term captures the effect ofself gravitation in the molecular cloud and ??2u accounts for the viscosityin the fluid.The Poisson equation,?2? = 4piG?,where G is the gravitation constant. The Poisson equation describes therelationship between the gravitational potential field and the distributionof mass density in the molecular cloud. The Poisson equation adds selfgravitation to the molecular cloud.The internal energy equation,??&?t+ P? ? (u) = 0,where & is the fluid?s internal energy per unit mass. The internal energyequation relates the pressure scalar field and the bulk velocity vector to theinternal energy. The internal energy equation is derived from the first lawof thermodynamics, which is an energy conservation statement. The zeroon the right hand side of the equation is a consequence of the cloud collapsebeing isothermal ( i.e., in an isothermal collapse Fourier?s law is equal tozero, because there are no temperature gradients).The induction equation,?B?t= ?? (u?B).The induction equation is a statement about the coupling of the chargedparticles motion in the cloud and the background magnetic field.The magnetic divergence constraint is? ?B = 0.102.1. The Physics of Star FormationThe magnetic divergence constraint physically states that magnetic monopolesare not treated.Lastly, the equation of state,P = c2s?,where cs is the sound speed. The equation of state describes how the pressureof the fluid relates to the density of the fluid, under isothermal conditions.Under adiabatic conditions the equation of state takes on the following formP ? ?? , where ? is the adiabatic index (i.e., a measure of the fluid?s degreesof freedom). The adiabatic sound speed is related to the adiabatic equationof state through the following relation.c2sa =?P??, (2.1)where csa is the adiabatic sound speed. Next, it is important to mention theassumptions used to construct the MHD theory. The implicit assumptionsthat one should be aware of when using the MHD model are as follows:? The gas components are not far from local thermodynamic equilib-rium.? Frequent collisions in the plasma keep the distribution functions, inlocal regions, approximately Maxwellian. ( First and second assump-tions are related and ensure that the system of particles can be treatedas a continuum ).? Small plasma parameter limit. Plasma parameter can be interpreted asa ratio of potential to kinetic gas energy. When the plasma parameteris small, the potential energy of a particle with the particles in itsneighbourhood is weak, in comparison to its kinetic energy, allowingfor free ions and electrons.? Low frequency Ampere?s law is assumed (i.e., displacement currentsare negligible), this assumption is reasonable when dealing with goodconductors, such as plasmas.? Time scales are large when compared to the cyclotron period (if thisassumptions does not hold then fluid is not coupled to the magneticfield and magnetic field can be ignored).The above assumptions may seem restrictive, but the MHD theory out-lined gives fairly reliable agreement with experiments, and leads to an un-derstanding of the essential fundamental physics involved. For a deeper112.2. Enzodiscussion of the MHD theory, and the assumption associated with it, seeref. [7].2.2 EnzoThe equations presented in Section 2.1 are a compact representationof the physics that unfolds in the extreme environments set in molecularclouds. When applying the MHD equations it is of interest to find out whatinitial conditions lead to (in)stabilities. To simplify the process of stabilityanalysis, linear stability theory is used. Linear stability analysis involves thefollowing algorithm:? Start with a stable solution (i.e., laminar flow) to the MHD equations.? Perturb the stable solution with a small, theoretically infinitesimal,perturbation.? Plug the perturbed solution into the governing physical equations andlinearize the equations by dropping all non-linear terms in the pertur-bation. At this point a first order Taylor approximation is used torelate perturbed physical quantities to initial physical quantities.? Solve the eigenvalue problem resulting from linearization of the per-turbation equations. Solutions that grow in space, or time are referredto as instabilities. Finally, solutions that are damped out are stableperturbations.Linear stability analysis offers important information about the condi-tions in which flows will change from laminar (i.e., well behaved flows)to destabilized flows that can lead to complicated turbulent behaviour. Byignoring the nonlinear terms, which describe the interactions of these per-turbations, linear analysis is restricted to the analysis of infinitesimal per-turbations. Therefore in a linearly unstable regime, linear theory can notpredict the outcome of the flow past a point where the perturbations growto a finite amplitude. For example, the Reynolds number of a flow gives ameasure of how stable the flow is to perturbations. Experiments and lin-ear stability analysis arguments show that as the Reynolds number of theflow increases, the flow changes from a predictable laminar flow, to unpre-dictable chaotic turbulent flow. When the Reynolds number is high enough,any solutions provided by linear analysis are not applicable, because the122.2. Enzonon-linear terms, initially ignored, become important. At this point infor-mation from experiments or numerical simulations must be used to ensureproper understanding of the flow?s nature.Star formation process involves many complex and unstable flows. Thismeans linear analysis can only describe star formation processes before per-turbations in the flow become finite in amplitude. Computational fluiddynamics (CFD) offers methods to study complex non-linear hydrodynam-ical processes, such as star formation. CFD techniques fall into two generalcategories Multi-grid, and classical numerical techniques. Multi-grid meth-ods provide a more robust and equivalently efficient counterpart to classicalnumerical techniques. Classical numerical techniques mainly fall into twolarge categories: direct methods and iterative methods. Direct methodsare highly precise, and efficient. However, they are specialized techniquesthat usually only apply to models arising from separable self-adjoint bound-ary value problems. On the other hand, iterative methods (or relaxationschemes) such as the Jacobi and Gauss-Seidel iterations, are easy to imple-ment and may be successfully applied to a more general set of problems thanmost direct methods.Single grid relaxation schemes have limitations; one important limita-tion is that these methods perform poorly at simultaneously dampeningboth high and low frequencies of the numerical error. Multi-grid methodswere constructed to address this limitation. Multi-grid methods use a coarsegrid solution to dampen the low frequency components of the error. Next,the coarse grid solution is used as an initial guess for the finer grid, which inturn dampens the high frequency error components. The use of multiple res-olution grids with iterative methods is competitive with fast direct methodsand is general, allowing for large set of applications. The following exampledemonstrates the benefits of using multiple grids. The problem solved inthe example is calculating the electric potential within an infinite squarecoaxial cable. Since there is no charge between the outer conductor and in-ner conductor, I use the Laplace equation to calculate the electric potential.For the purpose of this example I will assume that the outer conductor isset at a fixed potential of 2 V , and the inner conductor is set at an electricpotential of 100 V . The partial differential equation I have to solve is thetwo dimensional Laplace equation, which is expressed as?2?e?x2+?2?e?y2= 0,where ?e, x, y are the electric potential, and the two spatial coordinates,respectively. The relaxation method I use to solve this problem is based on132.2. Enzoa property of the Laplace equation, which states that any interior point in asolution of the Laplace equation must be the average of the points in a circlecentred on that same interior point. I exploit this property and iterated overthe grid multiple times replacing each interior point with the average of itsfour nearest neighbours. This was done iteratively until the change betweenthe previous grid iteration and the current solution was smaller then somepredetermined error threshold. To demonstrate the merits of multiple-gridmethods, I solved the problem in two different ways: single grid iterativemethod, and multiple-grid ( i.e., I used two grids for simplicity ) iterativemethod. For the single grid calculation the interior of the grid was initializedto the average of the boundaries. Similarly, the coarse grid for the multi-grid scheme was initialized in the same manner and its solution was fed intothe fine grid as an initial guess. The final resolution of both methods was1282, and the coarse grid had a resolution of 642. I ignored the numberof iteration it took to compute the coarse resolution solution, because thecomputational time was negligible when compared to the finer resolutioncomputational time. Figure 2.1 plots the two dimensional solutions for thepotential where the different colours represent different magnitudes of thepotential. The important thing to note from Figure 2.1 is that the twolast solution have a similar amount of error in them, but the multiple-gridmethod required fewer iterations. This highlights the main advantage of themultiple-grid method it generally converges faster to the solution then itssingle-grid counter part.Physical problems encountered in the field of star formation have largevariation in the order of magnitude of their spatial and temporal scale. Nu-merical techniques applied to these problems require highly refined spatialgrids to resolve accurate solutions. Often there are regions within the prob-lem?s domain in which the grid?s high resolution is not needed, and a coarsergrid would suffice. It is not efficient to allow local accuracy requirements dic-tate the global grid spacing. In these sorts of problems, highly refined meshare globally computationally wasteful, and can lead to computationally non-tractable problems with certain accuracy requirements. By locally refiningthe grid where it is needed, Adaptive Mesh Refinement (i.e., AMR) allowsoptimal spatial distribution of computational resources, while not compro-mising on desired accuracy. For example, the collapse of a molecular cloudpresented in Chapter 1 dealt with large variations in the Jeans length overthe molecular cloud?s domain. One technique employed by AMR is refine-ment of spatial resolution dependant on the Jeans length. Hence, the MHDequations remain locally linear, and finite difference methods are still usedwith high accuracy. Adaptive Mesh refinement employs multi-grid tech-142.2. EnzoFigure 2.1: The top three plots are of the Multi-Grid (MG) method forthree different iterations. The bottom three plots are of the Single-Grid(SG) method for three different iterations. The last plot for each respectivemethod has roughly the same measure of error, but the MG method hastaken fewer iterations.152.3. Post-Processingnique in which grids with different resolutions are used in conjunction withthe spatial refinement criteria, like the Jeans length. AMR is an extensionof the multi-grid techniques. The most important difference being that inAMR the finer grids do not cover the entire domain of the problem.Enzo is a 3D MPI-parallel Eulerian code block structured AMR code,used in my thesis to model the process of star formation. MPI-parallel refersto the message passing interface used by Enzo to allow for its algorithm tobe ran on large number of processors simultaneously. Eulerian means thatthe code computes all the physical quantities on a grid instead of movingwith the fluid. The code is implemented in C++ and Fortran 77. High-levelfunctions and data structures are implemented in C++ and computation-ally intensive lower-level functions are implemented in Fortran. Enzo usesthe HDF57 data format to write out data and restart files in a platform-independent format. Reference [17] introduces Enzo and gives an in-depthdiscussion of the numerical techniques used in its implementation.Enzo is complicated Physics code. It has numerous test problems, whichset the code up for different physical scenarios. In addition, each individualproblem has a plethora of parameters which control the particular physicalscenario?s initial configuration. The one specific problem my research wascentred around was called Hydrodynamics and MHD Turbulence problemsin Star Formation, and in Enzo this problem was coded as problem type 106.In order to run problem type 106 in Enzo, a problem specific parameterfile must be setup. There are a large variety of parameters that can beadjusted in the Turbulence in Star Formation problem. In our research weonly adjusted a subset of the parameters available, and used default settingprovided by the problem type for the rest of the inputs. For a detaileddescription of the inputs explicitly set in Enzo, refer to chapter 3.2.3 Post-ProcessingObservational data suited to the study of turbulent molecular clouds,is represented in the form of PPV cubes (i.e., position-position-velocitydata). Inherently, PPV cubes are a consequence of data recorded from radiotelescopes. Radio telescopes take 2D grid images of spectra in the sky. Eachspectrum is a plot of intensity versus frequency. The frequency information isused with the Doppler effect to convert the data from frequency to velocitymeasure. In this study we use CO line emission data. The CO emissionline is a result of rotational transition, corresponding to the molecule goingfrom a high rotation state to a lower one and emitting a photon in the162.3. Post-Processingprocess. Studies have shown that the abundance of carbon monoxide (CO)is sensitive to properties of the gaseous medium, such as density and velocity.For example, ref. [24] concludes that CO formation strongly depends onthe molecular cloud properties. These observational PPV cubes containinformation on the projected two-dimensional distribution of the gas andon the radial component of velocities. It is important to understand thatobservational PPV cubes are a result of complex physical process involvingintegration of the radiation intensity along the line of sight (LOS) through ahighly heterogenous density field. The complex projection process of turning3D data into PPV cube is discussed in ref. [21].Equations 2.2 and 2.3 describe how the frequency channels are convertedinto velocity channels.fc =f01 + z. (2.2)vc = cf20 ? f2cf20 + f2c, (2.3)where z, f0, fc, vc, v1, v2, and c are redshift, rest frequency of the emissionline (EL), central frequency of the EL, LOS velocity, velocities associatedwith both ends of the line profile, and speed of light, respectively.Observational data is represented in PPV cubes, but simulation outputis represented in PPP cubes (i.e., three dimensional space of fluid physicalproperties such as density, velocity, etc) . Therefore, simulation data mustbe converted to synthetic PPV cubes, for purpose of comparison.To compute the intensity of the CO line for the PPV cubes in the simu-lated turbulent molecular cloud, I use the radiative transfer code RADMC-3D. RADMC-3D calculates the population levels through the Large VelocityGradient approximation (LVG), which uses local velocity gradients to esti-mate photon escape probabilities in each zone of the simulation. The LVGapproximation was derived in ref. [25], and the approximation was extendedin ref. [18] to molecular clouds. The LVG implementation in RADMC-3Dis described in ref. [24]. In molecular clouds the assumptions of LVG ap-proximations do not hold up, but with the extension derived in ref. [18],and implemented in RADMC-3D, it has been shown that the extended LVGapproximation can be used in the context of molecular clouds, refer to [18]for more details.172.4. Comparison Statistic2.4 Comparison StatisticNumerical simulations are a chief tool used to investigate the effects ofnonlinear interaction. In the process of star formation, nonlinear effects areimportant and drive complex processes such as turbulence. By finding agood match between the observations and simulations, we recognize a setof possible physical parameters that drive the physical conditions withinmolecular clouds. The continual improvement in quality of simulations andobservations, gives rise to a demand for a framework by which comparisonstools for observations and simulations can be evaluated. Many statisticaltools are currently developed ad hoc, and it is unclear what each individualtool is telling us about the underlying physics. This thesis presents anelementary method in which the test statistics can be evaluated and theirphysical sensitivity probed ref. [28]. We assume that the statistic used is ascalar function belonging to the non-negative real numbers.A good distance metric for PPV cubes should satisfy several properties.The following list details important aspects that desirable metric distancestatistics should satisfy. Let O1and O2 represent two PPV data sets, andd(O1, O2) represent the distance metric function. The following list is aconcise description of the properties of a well behaved PPV metric, for amore detailed exposition refer to [28].? The distance metric should be invariant with respect to shifts in thecoordinates of the data. For example, d(O1, O2) = 0 if O1(x, y, v) =O2(x + ?x, y, v)? The statistic should not be sensitive to spatial ordering of pixels. Anystatistic that ignores the spatial relationships between pixels will fulfillthis criterion [28].? The statistic should be sensitive to spatial scaling (i.e., d(O1, O2) > 0if O1(x, y, v) = O2(sx, sy, v)).? Independent of observational noise.? Equality between different realizations of the same physical process.Practically, this means that d(O1, O2) ? 0 for simulation data withsame physical parameter setting, but different random seed (randomseed is used for driving the turbulence stochastically).We chose the following four statistics to demonstrate the Astro-statisticcomparison framework we have implemented. Our framework allows for182.4. Comparison Statisticsorting out what physical parameters each statistic is most sensitive to,and offers a consistent way in which multiple statistics are compared. Thefollowing subsections briefly explain the mechanics of the four statistics.2.4.1 Principal Component AnalysisPrincipal component analysis (PCA) is a multivariate statistical tool thatre-parameterizes the original explanatory variables into subset of orthogo-nal vectors. The PCA technique transforms the original variables in sucha way that the first principal component encompasses the most responsevariation; next, the succeeding variable is chosen such that it describes thesecond highest variation in the response while still being orthogonal to firstprincipal component. The same criteria applies to the rest of the principalcomponents.[11] first proposed applying PCA to spectral line data cubes of molec-ular clouds. Subsequently, PCA was shown to be sensitive to the structurefunction of the underlying turbulent flow [5]. PCA considers a PPV datacube as an ensemble of Nv objects each with Nxy attributes of intensity ateach velocity channel [11]. Each data PPV cube is turned into a matrixwith elements Oij , where the index i runs over all the spatial positions in aset velocity slice, and j is index varying between velocity slices. Nxy and Nvare the total of different spatial intensity pixels and number of velocity chan-nels, respectively [4]. Next, a covariance matrix, comparing the variancebetween different velocity planes, is constructed for both cubes. Equation2.4 describes the covariance matrix calculation.Sjk =1NxyNv?iOijOik. (2.4)Thereafter, we compute eigenvalues ? and eigenvectors u from Su = ?u,for S1 and S2 (where subscripts 1 and 2 correspond to the two data cubesO1 and O2). To construct a distance measure, we sort the two vectors ofeigenvalues in order of absolute value and use equation 2.5 to quantify thedistance.dPCA(O1, O2) =[?i(?1i ? ?2i)2?i ?1i?i ?2i]12, (2.5)where the sum over index i runs over the vector constructed from allthe eigenvalues found during the Su = ?u calculation, if the length of ?eigenvalue vector exceeds 50 elements then it is truncated.192.4. Comparison Statistic2.4.2 The Spectral Correlation FunctionThe Spectral Correlation Function (SCF) was first proposed by ref. [23].Subsequently, ref. [20] have shown that SCF is useful tool for detectingdifferences in turbulence properties when comparing simulation to observa-tional data. Different varieties of the SCF algorithm can be found in theliterature, we use SCF as a normalized root-mean-square difference betweenspectra separated by an offset vector l ref. [28]:S(l) = 1???? |O(r, ?)?O(r+ l, ?)|2?? |O(r, ?)|2 +?? |O(r+ l, ?)|2 , (2.6)where O and r are the PPV datacube and the vector representing the twospatial coordinates within the datacube, respectively.Equation 2.6 generates a surface for each PPV cube, and the two surfacesare compared using the metric displayed in equation 2.7.dSCF(O1, O2) =??l[S1(l)? S2(l)]2 (2.7)2.4.3 Cramer statisticThe Cramer statistic is a nonparametric two-sample-test that can beapplied on multivariate data as well as univariate data. The test statistic isthe difference of the sum of all the Euclidean inter-point distances betweenthe random variables from the two different samples, and one-half of the twocorresponding sums of distances of the variables within the same sample [2].Our implementation of the Cramer statistics calculates the statistic for twobivariate distributions. Restricting computation to bivariate distributionsallows for faster computational execution, when dealing with high resolutionmolecular cloud data. We flattened the spatial dimensions of the PPV cubesto compute the Cramer metric. To achieve flattening of the cubes, each twodimensional velocity slice is flattened into a vector and sorted by intensity.The bottom 75 precent of each velocity slice is discarded; this makes thestatistic more robust to the instrumental noise convolved into the observa-tional PPV data cubes. In the case of two data cubes with different spatialdimensions velocity slice intensities are truncated to the smaller cube. Thesimilarity value returned by Cramer statistic is dependent on the Intensityvalues in the PPV cubes; therefore, we normalize the intensities values in theflatten cubes (i.e., position-velocity matrices) by the spectral norm, which202.5. Experimental Designfor a real matrix is ||O||2 =?maxeigenvalues(OTO). We have found that aboveprescribed normalization keeps the value of the Cramer metric independentof intensity ranges. Equation 2.8 is the mathematical representation of theCramer metric.dCramer(O1, O2) =Nv1Nv2Nv1 +Nv2??2Nv1Nv2Nv1?i=0Nv2?j=0?(|| *O1i ? *O2j ||2)?1N2v1Nv1?i=0Nv1?j=1?(|| *O1i ? *O1j ||2)?1N2v2Nv2?i=0Nv2?j=1?(|| *O2i ? *O2j ||2)?? , (2.8)where Nv1 and Nv2 are the total number of velocity channels for O1 andO2, respectively. The ? function is the kernel function, our implementationof Cramer statistics includes several different kernel functions. In my thesiswe used the following kernel:?(z) =?z2. (2.9)The above kernel out performed the rest of the available kernels, when sub-jected to experimental testing. Lastly, it should be noted that our implemen-tation of Cramer statistic is not sensitive to pixel ordering, but does seemto satisfy the rest of the distance metric criteria, and produced interestingresults discussed in chapter 4.2.4.4 Maximum Mean DiscrepancyMaximum Mean Discrepancy (MMD) is the supremum of the differencein expectations (i.e., largest difference in expectations) over functions ina unit ball of a reproducing kernel Hilbert space (RKHS) [10]. The PPVsimulation data for the MMD statistic was prepared in the same way as theCramer statistic, described in Section 2.4.3. We used an MMD implemen-tation provided in kernlab, written in the statistical language R. We treatedMMD as a black box statistic, refer to [10] for a detail theoretical discussionof the MMD metric measure.2.5 Experimental DesignOne factor at a time (OFAT) experiment is an experiment where one fac-tor is varied while the remaining factors are held constant. OFAT designs are212.6. Response Analysisused regularly in the field in Astronomy, and more specifically in the study ofStar Formation. The following are several examples in which OFAT designswhere used in the literature, [24], [22], [1], [6]. An alternative to OFATexperiments is designed experiments, which vary multiple factors simulta-neously. There are many advantages to using designed experiments over thefield standard of OFAT experiments. The following list highlights the majordisadvantages of using OFAT designs as oppose to designed experiments.? Extra experimental runs are necessary when using designed experi-ments, but the estimates of the effects are more precise (i.e., lessvariable)? When dealing with experiments that investigate two or more factors,the estimates of interactions are not accessible; designed experimentsallow estimation of interaction effects.? Designed experiments intentionally probe a larger region of factorspace, which improves the prediction of the response surface in ques-tion.In this thesis we consider using two-level orthogonal factorial designs.For a deep discussion, refer to chapter 4.2.6 Response AnalysisThe response analysis is a major part of the research results presented.Broadly speaking, several techniques falling under the topic of multiple re-gression analysis were used, with variation on the underlying model errorassumptions. In addition, regression for both replicated and unreplicatedexperiments were employed. Replicates were not true replicates. Instead ofrepeating experimental runs, we used time steps from one simulation run asits pseudo-replicates. For replicated analysis the following techniques wereused ordinary least squares, generalized least squares, and weighted leastsquares. Different transformation of the response and factors were exper-imented with as well. In contrast, in the unreplicated experiment setting,three methods were used: one graphical, and two quantitative. The graph-ical technique was half normal plot of the coefficient estimates. Next, thetwo quantitative methods were, design centre point replication, and Lenth?smethod. For a lengthy discussion of the regression techniques used refer to4.22Chapter 3Methodology3.1 Physical Parameters/Factors ran in thesimulationThis section describes the physical factors varied in the simulations. Idevelop a physical intuition for these factors in order to better understandthe results presented in Chapter 4.3.1.1 Mach NumberSupersonic turbulence is fundamental to understanding star formationin molecular clouds. Highly compressible turbulence is believed to be aimportant driver of star formation in dense molecular clouds [12]. In oursimulations, we use the Mach number to vary the effects of turbulence on themolecular could. The Mach number is a unit-less ratio of the gas?s velocityover the sound speed in the cloud (i.e., isothermal sound speed is used).Higher Mach numbers M lead to stronger shocks and thus to higher densitycontrasts. For instance, the density jump in a non-magnetized, isothermalshock is proportional to M2 [9]. In our simulations, we use the MachNumber to effectively vary the effects of turbulence on the molecular could.The Mach number is a ratio of the gas?s velocity over the sound speed in thecloud. We used two settings for Mach number: 6 and 18. The crossing timeis defined as the time required for a gas to move from one side of the clusterto the other, if it is travelling at the typical random speed of objects in thecloud. This is the shortest time-scale on which significant hydrodynamicalevents can occur. In Enzo the crossing time parameter is called TimeUnits.Equation 3.1 describes how Mach number parameter is inputed into thesimulation.tU =LMcs, (3.1)233.1. Physical Parameters/Factors ran in the simulationwhere tU , L, M, and cs, are the crossing time, length of the simulatedmolecular cloud box (our simulation box sizes were 10 parsecs on each side),Mach number, and sound speed (refer to section 3.1.4 for details on cs),respectively.3.1.2 Magnetic FieldIn the setting of molecular clouds we are dealing with ionized plasma,which introduces magnetic field into the dynamics of the cloud. Studieshave shown that observations of molecular clouds are consistent with thepresence of super Alfvenic turbulence [12]. Magnetic tension is known tostabilize the hydrodynamics of unstable post-shock shear layers. Magneticpressure modifies shock jump conditions and affects final simulation starmass distribution when compared to predictions made based on cases withpure turbulence fragmentation [12]. In our simulations, we varied the Mag-netic field between two setting: 10?6 G, and 10?4 G.3.1.3 Driving scaleThere are three important spatial scales when discussing turbulence:driving scale, inertial scale, and dissipation scale. Driving scale is the size ofeddies for which energy is injected into the turbulent system. The eddies inthe driving scale obtain energy from the bulk flow of the fluid. For example,when one stirs their cup of coffee in the morning the initial size of the eddiesgenerated by the spoon are in the driving scale regime (i.e., energy injectionscale). The inertial scale is a subrange in which turbulent energy is passeddown from the driving scale down to the micro- scale, or dissipative scale.The inertial scale is the range in which power laws are stated. Power lawsare equations which describe the distribution of energy across the differenteddies, or scales. For example, a famous power law is Kolmogrov?s powerlaw.E(k)k2 ? &23 k?53 ,where E, &, k, are the energy distribution, the energy per unit mass per unittime transferred from the large eddies to small eddies, and k is the wavenumber, respectively. Lastly, the dissipative scale is the regime in whichviscosity become significant and dissipates the fluid?s kinetic energy. Thefollowing Lewis F. Richardson?s quote describes the three turbulent spatialregimes in a concise way:243.2. Physical Fixed Parameters ran in the simulationOde to turbulent Flow: Big whirls have little whirls that feed ontheir velocity, and little whirls have lesser whirls, and so on toviscosity.In our simulations the driving scale parameter adjusts the range corre-sponding to the driving scale. More specifically, we used an Enzo parametercalled CloudType that adjusts the maximum and minimum wave numberson which the turbulence energy injection is driven.3.1.4 TemperatureIn our simulations, we assumed that the molecular clouds are isother-mal. This assumption is reasonable, because the time scales on which equi-librium is reached are shorter than the hydrodynamical time scales. Thedifference in time scales allows the temperature of the cloud to adjust backto equilibrium before new hydrodynamical calculations are performed. Thetemperature parameter, adjusted in our simulation, appears in the Enzo pa-rameter list through the isothermal sound speed. For our particular problemtype in Enzo the sound speed has two parameter names which have to beset, EOSSoundSpeed and SoundVelocity. Equation 3.2 is simple equationdemonstrating how the temperature is used to calculate the isothermal soundspeed.cs =?kT?mh, (3.2)where k, ?, mh are mean mass of the gas relative to hydrogen particle,Boltzmann constant, and mass of hydrogen in grams, respectively.3.2 Physical Fixed Parameters ran in thesimulationAs mentioned in chapter 2, the Enzo parameter files have a large numberof possible inputs. For our particular simulation runs we used the defaultEnzo parameter file provided for problem type 106, and we varied inputsbased on factor setting described previously, in Section 3.1. There are severalparameter settings that we set specifically for our simulations, but were keptfixed across simulation runs. Table 3.1 highlights a set of important fixedparameters, inputed statically across all simulations. In post processingstage, we used 13CO as a tracer molecule. The statistics were computed253.3. Data Collection and Analysis Summaryusing the last 9 time steps of each simulation run. Time steps were recordedevery tenth of crossing time unit. Lastly, at this point we did not add noiseto the simulated PPV cubes. Noise will be added in future work to emulatethe instrumental noise present in observational data.Table 3.1: The parameters listed in this table were held fixed across simu-lation runs.Fixed Parameters Value DescriptionStopTime (Turbulence) 2 stop time in crossing units,turbulence mixing stage happen over the range 0-2StopTime (Gravity) 3 stop time in crossing units, turbulence andself gravity stage happen over the range 2-3TopGridDimensions 128 Grid resolution in simulated box(Left/Right)FaceBoundaryConditions 3 the boundary conditions for all cubefaces, 3 is coded for periodic boundary conditionsDensity 100 initial density in # of hydrogen atoms per ccLengthUnits 10 simulation box size in parsecsStarMakerMinimumMass 1 the minimum mass of star particle in solar massesMaximumRefinementLevel 3 # number of AMR refinement levelsStarParticleCreation 16 Sink particles: pure sink particle,or star particle with wind feedbackStaticHierarchy 0/1 set to 1 (static) for turbulence mixing stageand 0 (adaptive) when self gravity is turned on3.3 Data Collection and Analysis SummaryWe ran 16 simulations of ENZO version 1.9 (with modifications imple-mented by Stella Offner), and base problem 106 (i.e., MHD Turbulenceproblems in Star Formation). We varied Mach number, magnetic field, driv-ing scale, and temperature, but left the variables in Table 3.1 constant. Thesimulations were designed to emulate a star forming cloud undergoing grav-itational collapse. We ran all the simulations for 2 crossing times, withgravity turned off. Next, we turned on gravity, and ran the simulations for263.3. Data Collection and Analysis Summaryone more crossing time. We used sink particles and adaptive mesh. Wepost processed the data using RADMC-3D through one face and created adata cube with a velocity resolution of 500 cms . We did not add white noiseto mimic the effects of the telescopes. Noise addition will be addressed infuture work; in this dissertation comparison analysis was only done betweensimulations. Therefore, noise addition was neglected. Then we comparedeach cube at a time step to cubes from other physical parameters using thefour comparison statistics, PCA, SCF, Cramer, and MMD. Finally, we useda battery of statistical tests to analyze the four response variables. The sta-tistical tests determined the significant subset of parameters and interactionsdriving the difference between simulation cubes.An example of a 2D slice of a PPP cube for different run setting anddiffer face post-processing is shown in Figures 1.3b, 1.3a, 1.3d, and 1.3c.Examples of PPV cube slices at the zero velocity channel, across differentsimulation runs, is shown in Figure 3.1. Lastly, an example of a table withCramer statistic compares over the 16 runs for a specific time step is shownin Table 3.2.Table 3.2: Cramer statistic comparing a fiducial run?s specific time step toall the rest of the 15 simulations. Zero is coded for the low setting of eachparameter, and one is conversely coded for the high setting.Mach Bfield Driv Temp Cramer1 0 0 0 0 0.002 1 1 0 0 1.313 1 1 0 1 2.834 1 1 1 0 1.265 1 1 1 1 2.866 1 0 0 0 2.087 1 0 0 1 2.968 1 0 1 0 1.889 1 0 1 1 2.9510 0 1 0 0 0.6411 0 1 0 1 0.1312 0 1 1 0 0.5613 0 1 1 1 0.2914 0 0 0 1 0.8315 0 0 1 0 0.4016 0 0 1 1 0.76273.3. Data Collection and Analysis SummaryFigure 3.1: Examples of PPV slices extracted from four different simulationruns. Panel (a) is the fiducial run which all comparisons are made back to.Fiducial run has Mach number of 6, magnetic field of 1 ?G, default widerdriving scale, and a temperature of 10 K. Panel (b) has Mach number setto 18 and all the rest of the parameters as seen in the fiducial run. Panel(c) has the smaller driving scale, which corresponds to the high setting ofdriving scale parameter. Lastly, panel (d) shows a slice for high magneticfield of 100 ?G and all the rest of the parameters are set at their respectivefiducial settings. This image was generated in ref. [28].283.4. Experimental Designs3.4 Experimental DesignsHydrodynamical simulations are computationally expensive. Enzo, thecode we are running, can take on the order of days to complete one simula-tion. Currently, we are running the hydrodynamical code at a 1283 resolu-tion. This resolution has been chosen for its relative speed to accuracy ratiofor testing PPV comparison statistic framework. In future work the resolu-tion of the simulation will be increased to as high as 5123 and maybe even10243, for increased physical accuracy. Simulations at these resolution couldtake up to weeks to complete. Next, the size of the current parameter spaceis quite large, and in future work we would like to add factors and investigatepotential curvature in the parameter space. The computational complexityof the hydrodynamical code and the shear size of the parameter space beinginvestigated make any attempt at a brute force exploration of the parameterspace highly infeasible with current resources. This is where experimentaldesign shines and allows for optimal resource management and maximal in-formation retrieval. Using the extensively developed theory of experimentaldesign, one can ensure that experimental plan chosen will guarantee easyanalysis and allow for informative conclusions.In this thesis, we greatly borrowed upon the theory of fractional fac-torial design. We used both baseline and orthogonal parametrization ofthe designed experiments. For examples and theoretical discussion of theexperimental methods used in this dissertation, refer to chapter 4.29Chapter 4ResultsAs described in section 2.4, we used four statistics, PCA, SCF, Cramer,and MMD to compare PPV cubes in our experiments. The response vari-able was computed by comparing all runs to a specified baseline run (i.e.,all factors are at their low settings), referred to as the fiducial run. The pa-rameters investigated were Mach number, magnetic field, driving scale, andtemperature. Standard linear model approach with least square estimationis used to analyze the data. Several experimental designs and parameteriza-tions were analyzed to investigate which parameters? effects are significantwith respect to the four statistics.4.1 Baseline DesignThe design of an experiment is the first and most important step in anyparameter study. There are many available designs that perform well underdifferent conditions. In particular, we looked at designs for four factors,Mach number, magnetic field, driving scale, and temperature. Each factorwas assigned two levels, where levels refer to the number of setting for eachfactor, refer to Table 4.1 for specific settings (note: 0 is coded for thelow setting and 1 for the high setting). Three of the factors, Mach number,Table 4.1: Baseline factor settingsMach Bfield (G) Driv Temp (K)0 6.00 10?6 0.00 10.001 18.00 10?4 7.00 30.00magnetic field, and temperature are continuous factors, whereas driving fieldis a categorial variable. In this thesis we consider a one factor at a time(OFAT) design under baseline parameterization, full factorial designs andlastly a fractional factorial design. The later two designs are consideredunder orthogonal parameterization in the following Section 4.2.304.1. Baseline DesignThe differences between data cubes using one of the four different dis-tance measures (PCA, SCF, Cramer and MMD) will be treated as the re-sponse and all analysis of the data will proceed using a linear model. Thesaturated model is given byyi = ?0 +k?i=1?ixi +k?1?i=1k?j=i+1?i?jxixj+k?2?i=1k?1?j=i+1k?l=j+1?i?j?lxixjxl ++?1?2???kx1x2 ? ? ? xk + &i, (4.1)where y represents the response, x represents the different factors, ?0 isthe model intercept, ?i are the main effects, and the higher order ?1?2???correspond to interaction effects of increasing order. Lastly, &i is the modelerror (note: we assumed &i ? N(0,?2)). For simplicity, we can write theabove model in matrix form asY = X? + ",where Y, X, ?, and " are the response vector, the design or model matrix,vector of linear effects, and the vector of errors, respectively. The leastsquares estimates of the model coefficients are,b =????????b0b1...bk????????= (XTX)?1XTY. (4.2)For the OFAT design, the design matrix isX =??????1 0 0 0 01 1 0 0 01 0 1 0 01 0 0 1 01 0 0 0 1??????. (4.3)For the OFAT design a bi is interpreted as the expected change in the re-sponse as the factor i is moved from the low to the high setting. This314.1. Baseline Designinterpretation of effect estimates is typical in astronomy, as is the the use ofthe OFAT design in the star formation literature, for examples see ref. [24],ref. [22], ref. [1], and ref. [6].Using the last nine time steps as replicates (pseudo-replicates) Table 4.2shows the analysis of the data using the PCA as the distance measure.Table 4.2: One factor at a time baseline design results for PCA statistic. Thefirst column gives the effect estimates and the last column gives the p-value,or measure of significance. The smaller the p-value the more significant thecorresponding effect estimate.Estimate Std. Error t value Pr(>|t|)(Intercept) 0.0000 0.0490 0.00 1.0000Mach 0.7284 0.0693 10.52 0.0000Bfield 0.5800 0.0693 8.37 0.0000Driv 0.2340 0.0693 3.38 0.0016Temp 0.4444 0.0693 6.42 0.0000Table 4.2 shows that all the main effects are statistically significant andgives corresponding estimates. The main drawback of the OFAT design isthat it does not allow for the estimation of interactions between the pa-rameters. In order to consider interactions it is necessary to run a differentdesign.A two-level full factorial design is found by considering all possible levelcombinations between the k factors giving 2k runs (i.e., 16 runs for ourexperiment) and allows for the estimation of all the main effects and theinteractions between these factors. Table 4.3 shows the results from thisanalysis and Figure 4.1 shows a select group of interaction plots for PCA,SCF, Cramer, and MMD, respectively. Figure 4.1 and Table 4.3 clearlyshows that the interactions between factors are not negligible, and must beconsidered in the statistical analysis.The full factorial design is uneconomical when simulations with higherresolution are run, or factors with more levels are investigated. Fractionalfactorial designs offer an alternative by only using a select number of runsfrom the full factorial design while still retaining most of the informationavailable in the full factorial design. Fractional factorials designs underbaseline parameterizations can be found in ref. [15], but they producebiased estimators of the main effects when interactions are significant. Thisis clear when comparing tables 4.3 and 4.4 for the PCA statistic. The main324.1. Baseline Design0.00.51.01.5MachCRAMER0 1 Temp100.30.50.70.9MachPCA0 1 Temp100.020.060.100.14MachSCF0 1 Bfield01?1.0?0.50.00.5MachMMD0 1 Bfield01Figure 4.1: Several factors that interact across the different distance metrics.If the two lines are not parallel then an interaction between the labeledparameters is present.334.1. Baseline DesignTable 4.3: PCA full factorial saturated model linear regression results.Estimate Std. Error t value Pr(>|t|)(Intercept) 0.0000 0.0402 0.00 1.0000Mach 0.7284 0.0569 12.81 0.0000Bfield 0.5800 0.0569 10.20 0.0000Driv 0.2340 0.0569 4.12 0.0001Temp 0.4444 0.0569 7.82 0.0000Mach:Bfield -0.7612 0.0804 -9.47 0.0000Mach:Driv -0.3275 0.0804 -4.07 0.0001Bfield:Driv -0.2535 0.0804 -3.15 0.0020Mach:Temp -0.2078 0.0804 -2.58 0.0109Bfield:Temp -0.8638 0.0804 -10.74 0.0000Driv:Temp -0.3702 0.0804 -4.60 0.0000Mach:Bfield:Driv 0.3447 0.1137 3.03 0.0029Mach:Bfield:Temp 1.0383 0.1137 9.13 0.0000Mach:Driv:Temp 0.5146 0.1137 4.53 0.0000Bfield:Driv:Temp 0.4444 0.1137 3.91 0.0002Mach:Bfield:Driv:Temp -0.6332 0.1608 -3.94 0.0001effects estimates are different under the two designs. A bias is observed inall of the statistics main effects when comparing the full factorial design tothe fractional factorial design under baseline parameterization.In the majority of the analysis above most of the effects are judged asbeing significant. This is a sign of oversensitivity and is usually caused by aviolation of the linear model assumptions. Figures 4.2a and 4.2b show theresidual versus run plot and a normal QQ-plot for PCA. If the assumptionsare met one would expect to see a uniform scatter of points around zero forFigure 4.2a, and the points on the QQ-plot should fall very close to theline for Figure 4.2b. Each of these plots clearly indicate that this is notthe case and the linear model assumptions are grossly violated. In order toimprove the fitted model, I considered transformations of the response andalso considered using Generalized Least squares with an autoregressive errorstructure. None of these measures resulted in an adequate fit to the data.All statistical methods I tried led to similar error assumption violations.The only way in which the replicated analysis could proceed is to actuallyreplicate the full suite of simulations, which is not possible. Fortunatelythere is a vast literature on analyzing unreplicated factorial design using anorthogonal parameterization.344.1. Baseline DesignTable 4.4: Fractional factorial for PCA statistic under baseline parameteri-zation regression results.Estimate Std. Error t value Pr(>|t|)(Intercept) 0.2608 0.0410 6.36 0.0000Mach 0.4123 0.0580 7.11 0.0000Bfield 0.1325 0.0410 3.23 0.0016Driv 0.0332 0.0410 0.81 0.4193Temp -0.0615 0.0410 -1.50 0.1364Mach:Bfield -0.2280 0.0580 -3.93 0.0001Bfield:Temp -0.0562 0.0580 -0.97 0.3343Mach:Temp 0.4103 0.0580 7.07 0.0000Before proceeding to unreplicated analysis, I decided to check if therewas second order curvature in the response surfaces. If the response surface,for any of the four statistics, is curved, the linear fit I am using would beinadequate in both replicated and unreplicated analysis. In order to testsecond order, or quadratic curvature, I ran four replicates at the centre ofthe full factorial design. I used a lack a fit test, which required the centrepoints and is implemented in R, to check if the response surface displayedsigns of curvature. The test works by setting up a Null hypothesis thatthere is no curvature, and an alternative hypothesis, assuming that there iscurvature. The test spits out a p-value, which gives the probability of gettingthe lack of fit measure observed, or a worse one, under the Null hypothesis.Therefore, the smaller the p-value the less likely the Null hypothesis is true,which means there is evidence for curvature. From Table 4.5, it is clearthat there is no concrete evidence that curvature is present. Aside from theSCF statistic, which might need closer look for the possibility of responsecurvature, the rest of the statistics exhibit a linear response surface over theinvestigated parameter space. This means that my choice of a linear model issufficient for the current data. This supports the use of unreplicated designanalysis.354.1. Baseline Design?2 ?1 0 1 2?0.20.00.20.4pcaTheoretical QuantilesQ?Q residual plot(a) If the assumption of error normalitywas met, one would expect to see thepoints in the QQ-plot fall very close tothe line.5 10 15?0.20.00.20.4runspca(b) If the assumption for constant vari-ance was met, one would expect to seea uniform scatter of points around zero.This is not the case in this plot.Table 4.5: Lack of Fit p-values for the four different statistics, or distancemetrics. Smaller p-values mean that there is more evidence for curvature.Note, p-values are probabilities. Therefore, they range between 0 and 1.The only metric that might exhibit signs of curvature is SCF.Stat Lack.of.Fit (p-value)PCA 0.74Cramer 0.95SCF 0.055MMD 0.89364.2. Orthogonal Design4.2 Orthogonal DesignIn order to switch from replicated to unreplicated analysis I used anorthogonal parameterization in which the fiducial level (0) is transformed to(?1) and the high level is kept at a coding of +1. The first five columns ofthe orthogonal full factorial model matrix are given asX =?????????????????????????????1 ?1 1 1 11 1 1 1 11 ?1 ?1 1 11 1 ?1 1 11 ?1 1 ?1 11 1 1 ?1 11 ?1 ?1 ?1 11 1 ?1 ?1 11 ?1 1 1 ?11 1 1 1 ?11 ?1 ?1 1 ?11 1 ?1 1 ?11 ?1 1 ?1 ?11 1 1 ?1 ?11 ?1 ?1 ?1 ?11 1 ?1 ?1 ?1?????????????????????????????(4.4)where the first column represents the intercept and the remaining 4 columnsshow the level setting for the parameters. For example the third row showsthat Enzo would be run with Mach = 6, Bfield = 1 ?G, Driv = [1, 2] andTemp = 30 K. Also notice that rows 15, 16, 14, 10 and 7, respectively, arethe 5 runs in the OFAT design.When using this parameterization, XTX = nI where I is the identitymatrix; hence, this coding or parameterization is referred to as orthogonaldesign. In this case, the covariance matrix is diagonal, and the estimatedeffects ?? are independent. Note, under previous baseline parameterizationthe covariance matrix was not a diagonal matrix and the effect estimateswere not independent. Under this framework each parameter estimate isgiven by ??i = 12 (y?+ ? y??) where y?+ and y?? are the means of the responsesat the high and low levels respectively. Interpreting this back to differencesin the original data set is complicated. The consequence of changing codingenables independent estimates of the effects, at the expense of a direct linkbetween changing physical parameters and the resulting effects in the OFAT0,1 coding. The orthogonal design parameters can be roughly interpreted374.2. Orthogonal Designas a measure of how the similarity metric changes as one changes from thehigh level to the low level.4.2.1 Fractional Factorial Orthogonal DesignThe full factorial design provides a more nuanced interpretation of theresponses as a function of changing initial conditions. However, running afull factorial design can require substantial computational resources whenlarger studies are conducted, as described in Section 4.1. We ran the two-level, full factorial design, because we were using 1283 resolution simulations.In addition, we choose a model where there are only four factors and eachfactor only has two levels. There are situations where one would like to runsimulations with higher resolution and more factors and a continuous rangeof levels. Running a full factorial on these sorts of requirements would beprohibitively expensive computationally. It is possible to run a subset ofthe full factorial design with the ability to estimate both the main effectsand some of the interactions. In many situations, prior knowledge can beused to identify a limited subset of the interactions that are likely to beimportant. Starting from the full factorial design of 16 runs, we select asubset of 8 runs that correspond to the fractional factorial design. Thisresults in the following 8 runs, which are commonly referred to as a 24?1fractional factorial design.X =????????????1 ?1 ?1 ?1 ?11 1 ?1 ?1 11 ?1 1 ?1 11 1 1 ?1 ?11 ?1 ?1 1 11 1 ?1 1 ?11 ?1 1 1 ?11 1 1 1 1????????????. (4.5)Again the columns 2 through 5 are used to assign the level setting for eachof the four factors. Since we have reduced the run size by a factor of twowe lose the ability to estimate 8 of the 15 effects that were estimated usingthe full factorial design. Of the 7 remaining degrees of freedom, one canestimate the effects of one of two types of terms. The main effects or the asubset of 3 interactions that in our case are Mach:Bfield, Mach:Temp, andBfield:Temp, as was done previously using the baseline fractional factorialdesign. The difference of using the orthogonal fractional factorial as opposed384.3. Unreplicated Design Analysisto the baseline fractional factorial is that the orthogonal design does not havebias in the main effects estimates.4.3 Unreplicated Design Analysis4.3.1 Normal Plots/Half-Normal PlotsNormal and half normal probability plots use the same statistical con-cept to pick out significant effect estimates. The only difference is how thetwo methods present their results. Half normal plots have a slight visualadvantage, stemming from the arbitrariness of the estimate?s sign when sig-nificance is tested, further explanation follows. The normal probability plotis most commonly used for assessing the normality assumption ( i.e., theQQ-plot used in previous section). In the context of orthogonal unreplicatedexperiments, [8] proposed the use of normal plots as a test or visual judge ofthe significance of the factorial effects. Let ?1 ? ... ? ?L be the sorted valuesof the effect estimates of a specific model, where L is the number of effectestimates. We plot these estimates against their corresponding coordinateson the normal probability scale:?i against ??1([i? 0.5]L)fori = 1, ..., L, (4.6)where ? is the cumulative distribution function of a standard normal randomvariable. Next, we graph a line through the points with ?is roughly equalto zero. Finally any effect whose corresponding points fall off the line aredeemed significant.The justification of above algorithm is as follows. We first assume thatthe estimates are normally distributed with means equal to the effect esti-mates. The assumption of normality is not overly restrictive, because theeffect estimates are averages of n observation (note: in our case the fullfactorial has n = 16, and the fractional factorial has n = 8) and the centrallimit theorem ensures that averages are nearly normal for n as small as 8.The null hypothesis is that all the effect estimates are zero. Under the null,if all the estimated effects have the same mean of zero, then all the pointsshould fall roughly on a straight line. Therefore, the normal probabilitydistribution is testing whether all the estimated effects have the same dis-tribution and any points falling of the line suggest that their correspondingestimates are active with respect to the rest of the estimates. This clevermethod does not require ?, but does implicitly require ? to be constantacross the estimates. In addition, drawing the line through the non-active394.3. Unreplicated Design Analysispoints and choosing active effects opens the analysis to subjective step inthe absence of a quantitative measure of ?. The half normal probability plotfollows almost the exact same procedure, but instead of sorting the effectestimates it sorts the magnitude of the effect estimates. In addition, themagnitude of the estimates is plotted against the coordinates of the halfnormal probability scale:|?i| against ??1(0.5 + 0.5[i ? 0.5]L)fori = 1, ..., L. (4.7)The advantage of the half normal probability plot is that the active effectsare all above the line, which is just a visual convenience. For a more detailedexamination of these two visual methods refer to [27].I used half normal plots to pick out the effects that are most importantin both the fractional factorial and full factorial designs. Figures 4.2, and4.3 are plots of fractional factorial and full factorial models, respectively.One thing to note from figures 4.2, and 4.3 is that the top two effects inthe full and fractional factorial designs are consistent. The agreement ofthe two designs half normal plots re-assures the use of the more economicalfractional factorial design. Half normal plots give a rough idea of whicheffects are active, but it is hard to objectively pick out the exact subset ofimportant effects. There are two more quantitative methods I will present,that address the lack of objectivity in the half normal plots: centre pointestimation of the error, and Lenth?s method.4.3.2 Centre Point Variance EstimationCentre point estimation of the error uses replication of the design?s centrepoints to estimates variance, which can then be used to perform the usualanalysis of the model. Tables 4.6 and 4.7 give a summary of the significanteffects found using centre points analysis.Table 4.6: Significant effects in Fractional Factorial Design.PCA Cramer SCF MMDMach Mach Nothing MachMach:Temp Temp Mach:TempMach:TempPCA, Cramer and MMD statistics are mostly consistent across fractionaland full factorial designs, whereas some information is lost for SCF moving404.3. Unreplicated Design Analysis0.5 1.0 1.50.00.20.40.6Half?normal quantilescramer effectsMach:TempTemp0.2 0.4 0.6 0.8 1.0 1.2 1.40.000.100.20Half?normal quantilespca effectsMach:TempMach0.2 0.4 0.6 0.8 1.0 1.2 1.40.000.010.020.030.04Half?normal quantilesscf effectsMach:BfieldMach0.5 1.0 1.50.000.040.080.12Half?normal quantilesmmd effectsMach:Temp TempFigure 4.2: Effect estimates half normal plot for fractional factorial model.The dots that are sitting high above the rest or effectively outliers are con-sidered as significant effects. I have labeled the two most significant effects.414.3. Unreplicated Design Analysis0.0 0.5 1.0 1.50.00.20.40.6Half?normal quantilescramer effectsMach:TempTemp0.0 0.5 1.0 1.50.000.100.20Half?normal quantilespca effectsMach:TempMach0.0 0.5 1.0 1.50.000.010.020.030.04Half?normal quantilesscf effectsMach:BfieldMach0.0 0.5 1.0 1.50.000.040.080.12Half?normal quantilesmmd effectsMach:TempTempFigure 4.3: Effect estimates half normal plot for the full factorial model. Thedots that are sitting high above the rest or effectively outliers are consideredas significant effects. I have labeled the two most significant effects.424.3. Unreplicated Design AnalysisTable 4.7: Significant effects in Full Factorial Design.PCA Cramer SCF MMDMach Mach Mach MachMach:Temp Temp TempMach:Temp Mach:Tempfrom the full factorial to the fractional factorial design. Running a regressionusing centre point replication is a purely quantitative method of selectingsignificant effects, which requires costly centre point replication. Lenth?smethod offers an alternative objective technique to analyze the unreplicatedexperiment, without the need for centre points.4.3.3 Lenth?s MethodLenth?s method uses quantitative analysis to select active factor, butallows for the error to be estimated without the need for centre points. [14]considered a robust estimator of the error for the effect estimates. Lenthcoined the term pseudo-standard error which he referred to as PSE:PSE = 1.5 ?median|??i<2.5s0||??i|, (4.8)where the ??is are the effect estimates, and s0 is defined in equation 4.9.s0 = 1.5 ?median|??i|. (4.9)PSE is used to obtain effect significance. The ?i are divided by PSE tocompute a t-like statistics. Effects are declared significant if the absolutevalue of the t-like statistic exceeds a critical value. The critical values aresimulated values taken from ref. [27], and correspond to 95-precent confi-dence interval for the 16 and 8 run unreplicated designs.To display the Lenth?s method results I have used an intuitive graphicalmethod presented in Lenth?s original paper of 1989. Effect estimates aregraphed in a bar graph with confidence intervals drawn. Figures 4.4, 4.5, 4.6,and 4.7 present the Lenth?s method results.434.3. Unreplicated Design AnalysisEffect EstimatesPCA?1.0?0.50.00.51.0MachBfieldTempDrivMach:BfieldMach:TempBfield:TempEffect EstimatesSCF?0.06?0.020.020.06MachBfieldTempDrivMach:Bfield Mach:TempBfield:TempFigure 4.4: PCA and SCF Lenth results for fractional factorial design. Note,that any bar protruding out of the two horizontal lines is considered signifi-cant. The significant main effect for PCA is Mach and the significant maineffects for SCF are Mach and Temp. PCA has lost the temperature signif-icance moving between the full and fractional factorial, see Figure 4.6 tohighlight the difference.Tables 4.8 and 4.9 summarize the significant effects. The fractionalfactorial results based Lenth?s method are fairly consistent with the centrepoint analysis for full factorial design. It is interesting to note that the resultsproduced using Lenth?s method on the full factorial design suggests thatinteraction effect with magnetic field and Mach number is also important.Finally, after analyzing the non-replicated data, using centre points, itis apparent that no evidence was found to suggest that strong curvature ispresent in our current parameter space. The consistency apparent in theunreplicated analysis suggests that the use of the pseudo-replicates (PRs)should be avoided. The PRs variation underestimate the inherent variationin the simulations causing all the effects to show up as significant.444.3. Unreplicated Design AnalysisEffect EstimatesCramer?0.4?0.20.00.20.4MachBfieldTempDrivMach:BfieldMach:TempBfield:TempEffect EstimatesMMD?0.4?0.20.00.20.4MachBfieldTempDrivMach:BfieldMach:TempBfield:TempFigure 4.5: Cramer and MMD Lenth results for fractional factorial design.Significance is judged as described in previous caption. Note, the significantmain effect for MMD is Mach and the significant main effects for Cramerare Mach and Temp. MMD has lost the temperature significance movingbetween the full and fractional factorial, see Figure 4.7 to highlight thedifference.454.3. Unreplicated Design AnalysisPCA?0.15?0.10?0.050.000.050.100.15Effect EstimatesMachBfieldDrivTempMach:Bfield Mach:DrivBfield:DrivMach:TempBfield:TempDriv:Temp SCF?0.06?0.04?0.020.000.020.040.06Effect EstimatesMachBfieldDriv TempMach:BfieldMach:DrivBfield:DrivMach:TempBfield:TempDriv:TempFigure 4.6: PCA and SCF Lenth results for full factorial design. Note, thatany bar protruding out of the two horizontal lines is considered significant.The significant main effects for both PCA and SCF are Mach and Temp. Iomitted the higher order interactions to make the plot more readable; theywere not significant.464.3. Unreplicated Design AnalysisCramer?0.2?0.10.00.10.2Effect EstimatesMachBfieldDrivTempMach:BfieldMach:DrivBfield:DrivMach:TempBfield:TempDriv:TempMMD?0.050.000.05Effect EstimatesMachBfieldDrivTempMach:BfieldMach:DrivBfield:DrivMach:TempBfield:TempDriv:TempFigure 4.7: Cramer and MMD Lenth results for full factorial design. Signifi-cance is judged as described in previous caption. Note, the main effects thatare significant for Cramer and MMD are the same as for PCA and SCF. Iomitted the higher order interactions to make the plot more readable; theywere not significant.474.3. Unreplicated Design AnalysisTable 4.8: Active Factors in Fractional Factorial Design using Lenth?sMethod.PCA Cramer SCF MMDMach Mach Mach MachMach:Temp Temp TempMach:Temp Mach:TempTable 4.9: Active Factors in Full Factorial Design using Lenth?s Method.PCA Cramer SCF MMDMach Mach Mach MachTemp Temp Temp TempMach:Bfield Mach:Temp Mach:Bfield Mach:BfieldMach:Temp Mach:TempBfield:TempMach:Bfield:Temp48Chapter 5Conclusion/Discussion5.1 Result SummaryInitially, I collected data for a 2-level full factorial design under base-line parametrization. I set out to show that running a one factor at atime design ignores important two-factor interactions, resulting in lackingmodels, with inadequate predictive power. I assumed that measurementson one simulation run across time steps would provide adequate predictionof error across simulations (i.e., used the time steps in one simulation as?pseudo-replicates?); in addition, the use of a 2-level design imposes a first-order model assumption ( i.e., no curvature can not be estimated). It wasshown that interaction effects were in fact significant, refer to figure 4.1for graphical evidence. In addition I demonstrated that fractional factorialdesigns under baseline parameterizations resulted in biased estimates of pa-rameters? effects. At this point I introduced the orthogonal parametrizationand demonstrated its benefits: unbiased estimates from fractional factorialdesign, and the many analysis tools develop to treat such designs. Next, Ipresented post-hoc checks of the underlying assumptions for the regressionused to analyze the different designs under orthogonal parameterization. Allthe assumptions were violated indicating that one of the original assump-tions about the data were not reasonable. After running additional centrepoint simulation runs the curvature assumption was tested. Originally I ex-pected that the response surface would be non-linear, but after running lackof fit tests using the centre points no evidence was found to suggest localsecond order curvature. Lastly, the pseudo-replicates were averaged overto convert the data into an unreplicated design. Three methods were usedto analyze the unreplicated design: half normal plots, centre point repli-cation, and Lenth?s method. All three methods showed consistent resultsacross PCA, Cramer, MMD, where as SCF showed weaker agreement withthe other three statistics, but still displayed a similar trend. I suspect SCFmight have issues with second order curvature. The lack of fit p-value forSCF was roughly 0.05, indicating that SCF?s response surface might havenon negligible second order curvature.495.2. Result Interpretation5.2 Result InterpretationUnreplicated analysis of the four distance metrics consistently high-lighted three significant effects: Mach number, temperature and their in-teraction effect. The significant effects are intuitively consistent with thesimulations? input process, which is evidence that the statistics capture realsimulation physics. Temperature is not directly inputed into Enzo; it in-fluences the simulation through isothermal sound speed (ISS). In addition,ISS affects the crossing time which is another direct Enzo input. Crossingtime is affected by both Mach number and ISS. Equation 5.1 describes theproportionality relation for crossing time, dependent on Mach number andtemperature.tU ?1M?T. (5.1)From equation 5.1 it is noted that temperature and Mach number bothaffect crossing time inversely. Hence, Mach number and temperature shouldaffect the response statistics in the same direction. In addition, the interac-tion effect should have the same direction as Mach number and temperature.Finally, temperature should affect the response statistic more gradually thanMach number, because of the square root relationship it has with the ISS(i.e., the magnitude of the effect estimate for temperature should be smallerthan Mach number effect size). The results produced with the unreplicatedanalysis are consistent with above theoretical intuition, refer to figures 4.4-?? for demonstration. It is interesting to note that PCA, SCF, and MMDfull factorial experiments predict a significant negative effect associated withthe magnetic field Mach number interaction. I expected magnetic field toshow up as significant. I know that if I visually inspect our PPV cubes,where magnetic field is set high, that the cubes look significantly differentthan the fiducial run?s PPV cube. The lack of magnetic field significanceis interesting and I believe that it arises from the way the PPV cubes arecollapsed before they are inputed into the statistics. Magnetic field is a fac-tor that shows up visually through the spatial structure of the PPV cube,and our current data collapse methods seem to destroy the spatial structureinformation. Lack of spatial structure may be the cause of lack of magneticfield significance, which is an interesting result on its own. My research teamand I are in the process of further investigating the potential causes of lackof magnetic field significance.505.3. Future investigation of Curvature5.3 Future investigation of CurvatureIn chapter 4 I assumed a first order model, or linear model; therefore,I did not estimate second order terms relating to curvature. If curvatureis in fact present, then inference based on the model is incorrect. I alsopresented evidence to suggest that second order curvature is not presentin our parameter hyper-cube for PCA, MMD, Cramer. In contrast, SCFshowed signs of plausible second order curvature. Further investigation ofSCF parameter space is needed to conclude wether curvature is significant.For the purpose of this dissertation, an implicit assumption about curvaturewas made. I decided to only test the design for the existence of secondorder, or quadratic curvature. The pitfall of the specific lack of fit testused in this paper, is that it can not distinguish between a truly linearrelationship and relationship exhibiting higher orders of non-linearity (i.e.,higher than second order). Figure 5.1 gives an example of the aliasing thatoccurs between linear and higher order curvature effects. Curves labeled byone and three would be considered both linear under our current lack of fittest. It is also possible that in the cases of PCA, Cramer, and MMD, I arein a regime of local linearity. If our response surfaces are only locally linear,extrapolation predictions based on the linear model would produce poorresults. In future work, we will further probe curvature using two methods:sensitivity analysis, and a more flexible design such as a latin hyper-cube.Sensitivity analysis will allow us to test the robustness of linear model,shedding light on wether or not the assumption of linearity is appropriatefor the purpose of extrapolation. Finally, using a latin hyper-cube design,will allow us to estimate curvature, while still keep the number of simulationruns down to a minimum.5.4 Future Research DirectionsThere are several categories that our research can be extended to in fu-ture work. This dissertation?s main goal was to develop a proof of concept,outlining a framework in which a suite of statistics for comparing PPV datacubes could be analyzed. In this paper, I used simulations that were runat a relatively low resolutions. An interesting question to answer is, whatresults will be produced as the simulations? resolutions are increased. Lowresolution simulations introduce noise into the final simulation solution, im-pairing the quality of the physics captured. Next, it would be of interest toadd different molecular tracers to the post processing stage of our current515.4. Future Research Directions111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111110.0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.[binary].1: Example of higher order curvature. The test I ran for curva-ture, tests if the centre point falls on a line drawn through the edges ofthe parameter space. This is a good test for second order curvature, butif the surface actually exhibits third order curvature the test might fail atdetecting curvature. This can be better understood from the above plot.525.4. Future Research Directionsframe work. Different tracers capture different information ( e.g., differentdensity profiles) about the molecular cloud, and statistics built on multi-ple tracers would add robustness to the similarity metrics. In our recentwork, I are only analyzing information produced from one face on the threedimensional PPP simulation cube. A natural extension of the statistics,would integrate all three faces of the simulation PPP cubes. In chapter 4 Ireached the conclusion that the pseudo-replicates used in the initial analysisperformed poorly at estimating the underlying error in the simulation runs.I suspect that the pseudo-replicates under-estimated the variance, becausethe final time steps in the simulation are in a regime where self gravity is en-gaged. Pseudo-replicates used from the earlier stage of the simulation wherethe molecular cloud is turbulent and mixing without gravity, might providegood estimator of the simulation variance; hence, allowing for accurate in-ference on the parameters, without the need for additional computationallytaxing true replicates. Finally, the grand objective for future work would beto compare high resolution cubes to actual observational data. Comparingsimulation to observation would allow for calibration of inputs to hydro-dynamical codes and for a greater understanding of the underlying physicsbringing about the process of star formation.53Bibliography[1] Attwood, R. E., Goodwin, S. P., Stamatellos, D., and Whitworth, A. P.(2009), ?Simulating Star Formation in Molecular Cloud Cores. IV. TheRole of Turbulence and Thermodynamics,? The Astrophysical Journal.[2] Baringhaus, L. and Franz, C. (2004), ?On a New Multivariate Two-Sample Test,? Journal of Multivariate Analysis.[3] Bate, M. R. (2009), ?Stellar, Brown Dwarf and Multiple Star Prop-erties from Hydrodynamical Simulations of Star Cluster Formation,?MNRAS.[4] Bertram, E. (2012), ?Statistical Analysis of Turbulence in the ISM,?Master?s thesis at Heidelberg Institute for theoretical Astrophysics.[5] Brunt, C. and Heyer, M. (2002), ?Interstellar Turbulence. I. Retrievalof Velocity Field Statistics,? ApJ.[6] Burkhart, B., Falceta-Gonalves, D., Kowal, G., and Lazarian, A.(2009), ?Density Studies of MHD Interstellar Turbulence: StatisticalMoments, Correlations and Bispectrum,? The Astrophysical Journal.[7] Choudhuri, R. A. (1998), The Physics of Fluids and Plasmas an intro-duction for astrophysicists, Cambridge University Press.[8] Daniel, C. (1959), ?Use of Half-Normal Plots in Interpreting FactorialTwo-Level Experiments,? Technometrics.[9] Federrath, C. and Klessen, S. R. (2012), ?The Star Formation Rateof Turbulent Magnetized Clouds:Comparing Theory, Simulations, andObservations,? DRAFT.[10] Gretton, A., Borgwardt, M. K., Rasch, J. M., and Scholkopf, B. (2012),?A Kernel Two-Sample Test,? Journal of Machine Learning Research.[11] Heyer, M. and Schloerb, F. (1997), ?Application of Principal Com-ponent Analysis to Large Scale Spectral Line Imaging Studies of theInterstellar Medium,? ApJ.54Bibliography[12] Kritsuk, G. A., Norman, L. M., Padoan, P., and Wanger., R. (2007),?The Statistics of Supersonic Isothermal Turbulence,? ApJ.[13] Larson, R. B. (1981), ?Turbulence and Star Formation in MolecularClouds,? Monthly Notices of the Royal Astronomical Society.[14] Lenth, V. R. (1989), ?Quick and Easy Analysis of Unreplicated Facto-rials,? Technometrics.[15] Mukerjee, R. and Tang, B. (2011), ?On Fractions of Two-Level FactorialUnder a Baseline Parameterization,? BIOMETRIKA.[16] NASA (2013), ?NASA?s Spitzer Science Centre,?http://www.physicsoftheuniverse.com.[17] O?Shea, W. B., Bryan, G., Bordner, J., Norman, L. M., Abel, T., Hark-ness, R., and Kritsuk, A. (2004), ?Introducing Enzo, an AMR Cosmol-ogy Application,? ARXIV.[18] Ossenkopf, V. (1997), ?The Sobolev Approximation in MolecularClouds,? New Astronomy.[19] Padoan, P., Goodman, A. A., and Juvela, M. (2003), ?The SpectralCorrelation Function of Molecular Clouds: A Statistical Test for The-oretical Models,? The Astrophysical Journal.[20] Padoan, P., Rosolowsky, W. E., and Goodman, A. A. (2000), ?TheEffects of Noise and Sampling on the Spectral Correlation Function,?The Astrophysical Journal.[21] Pichardo, B. e. a. (2000), ?On The Effects Of Projection On Morphol-ogy,? The Astrophysical Journal.[22] Rosas-Guevara, Y., Vazquez-Semadeni, E., Gomez, C. G., and Jappsen,A.-K. (2010), ?Dependence of The Star Formation Efficiency on GlobalParameters of Molecular Clouds,? The Astrophysical Journal.[23] Rosolowsky, W. E., Goodman A. A., 2, . W. J. D., and J., W. P.(1999), ?The Spectral Correlation Function: A New Tool for AnalyzingSpectral Line Maps,? The Astrophysical Journal.[24] Shetty, R., Glover, C. S., Dullemond, P. C., and Klessen, S. R. (2011),?Modelling CO Emission - I. CO as a Column Density Tracer and TheX Factor in Molecular Clouds,? The Astrophysical Journal.55Bibliography[25] Sobolev, V. V. (1957), ?The Diffusion Radiation in Nebulae and StellarEnvelopes,? Soviet Astronomy, Vol. 1, p.678.[26] STScI and NASA (2013), ?Hubble Site,? http://hubblesite.org.[27] Wu, C. F. J. and Hamada, S. M. (2009), Experiments: Planning, Anal-ysis, and Optimization, Wiley.[28] Yeremi, M., Flynn, J. M., Loeppky, J., Rosolowsky, E., and Offner, S.(2013), ?Comparing Astrostatistics,? in Preperation.56
Thesis/Dissertation
eng
Interdisciplinary Studies
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Statistically calibrating comparison metrics of star formation simulations.
Text