Far-infrared (FIR) fine structure lines (FSLs) are important tools for studying galaxies nearby and far away, and observations have made progress in revealing star bursts hidden in nearby dusty galaxies, as well as confirming or discovering hundreds of galaxies in the early universe. But interpretation of FIR FSL observations remains a challenging task, especially when galaxy integrated gas emissions from both ionized and neutral gas are involved. This thesis aims to provide a comprehensive and systematic framework for understanding FIR FSLs in both observational and theoretical manners. In the first part, we present comprehensive catalogs of the global FIR FSL data of both low- and high-z galaxies collected from the literature. We calibrate the commonly used diagnostic relations on integrated data and explore new empirical relations among FIR FSL luminosities. We find a strong dependence of FSLs on elemental abundance, suggesting that FIR FSLs are primary tracers of the bulk mass of the emitting ion. We also find non-trivial active galactic nuclei (AGN) contribution to the [O III]88μm line emission. By combining spectral lines of various wavelengths and origins, we find very tight correlations between the neutral gas and ionized gas emission, as well as a concordance between the optical and FIR line ratios. This makes FIR FSLs excellent tracers for the total amount of the emitting ions. Furthermore, we find the [C II] "deficit" problem is just part of a fundamental problem, where all lines show the same "deficit" trend including neutral or ionized gas emissions, and even extinction corrected H⍺, while the heating, properties, and luminosity of dust decouples from that of gas, resulting in a gas-dust dichotomy. We argue that this could be caused by an infrared excess or a deficit of ionized gas emission. New pictures of distribution/origin of the line and/or dust emission must be developed to reconcile the observational data, and we discuss some promising scenarios. This problem also has profound implications on the applicability of the photodissociation region (PDR) model and the utility of various tracers of the star formation rate (SFR). In the second part of the thesis, we use photoionization models to study gas emission from both ionized and neutral gas in a coherent structure. We demonstrate that a power-law fit on model parameters can quantitatively describe the scaling of most FIR FSL ratios on model physical properties. The known line ratios are investigated with respect to density, radiation field strength, elemental abundance, and electron temperature. With the help of grid diagrams and power-law fits, we confirm the primary dependence of the line ratios on the relevant physical quantities, and reproduce most of the empirical relations. We discuss the general properties of the ISM inferred from line diagnostics, and argue that a strong parameter marginalization exists in real galaxies in terms of the density and the O/H-U-Q1/Q0 correlation. Challenges also emerge, including the degeneracy between the ionization parameter and the radiation hardness, neutral gas density in models, the difficulty in modelling a realistic dusty H II region, etc. In the third part of the thesis, we introduce a data reduction package developed to process the observational data of ZEUS-2, a FIR spectrometer built to study FIR FSLs in high redshift galaxies. New tests are presented to characterize instrument performance along with explanations based on a simplified transition edge sensor (TES) theory that focuses on steady states. The new noise removal processes based on independent component analysis (ICA) is demonstrated. We point out the current bottleneck limiting the instrument performance is likely the telescope efficiency. An updated estimation of instrument sensitivity is presented using the best-informed knowledge of measured noise and efficiency. In the end, we summarize the observational and theoretical limitations of FIR FSL studies to date, and discuss the prospects of these unique lines.