Variability of Blazars and Blazar Models over 38 Years

Since 1978, when BL Lac objects and violently variable quasars were married to become “blazars”, physical interpretations of the phenomenon have evolved. Remarkably, though, the general picture of relativistic jets beaming their radiation in our direction, proposed that year by Blandford and Rees, remains intact. The main stress on theoretical models has come from observations that reveal ever more extreme variability requiring rampant particle acceleration on time-scales of minutes, often parsecs away from the central black hole. Here the author reviews many of the observations and theoretical ideas that have shaped his studies of blazars over about 40 years. This leads to his preferred scenario that blazar jets contain a helical magnetic field close to the black hole, turbulent plasma on parsecs scales, and both standing and moving shock waves. Particle acceleration can then occur in multiple stages involving the second-order Fermi process, magnetic reconnections, and modest jumps in energy at shock fronts. The most extreme variability, as well as brightness temperatures ~100 times the inverse Compton limit, probably require occasional exceptionally high bulk Lorentz factors. These can result, for example, from supersonic, relativistic turbulence, or ultra-relativistic flows propelled from sites of magnetic reconnection. Future efforts in these and other areas can determine whether these potential solutions are valid.


Introduction
It has now been 38 years since the study of rapidly variable extragalactic objects received the focus provided by a name: "blazars", coined by Ed Spiegel during his post-banquet monologue at the Pittsburgh Conference on BL Lac Objects [1].Within a few years, initially reluctant journal editors relented and allowed the terminology to appear in the literature when referring to BL Lacertae objects and violently variable quasars.Most of us now consider any radio-bright active galactic nucleus with a relativistic jet pointing within about 10 • of the line of sight to be a blazar.The author has been involved in attempts to explain blazars for about 40 years, during which time the observed variations in flux and polarization have become even more extreme.As is often the case, nature has proven to be much more exotic than anticipated, leaving us to scramble to interpret events that seem at first to defy its laws.Yet, ironically, most of our current explanations are still based on the relativistic jet model proposed by Blandford and Rees at the conference in 1978 [2].After almost 40 years, our understanding of the physical processes that occur inside such jets remains incomplete.
Here the author provides a retrospective of the evolution of blazar research over his career, emphasizing events and ideas that have shaped his thinking.Because of page limitations, this is a selective review that omits many important contributions to the field, for which he begs forgiveness.He hopes that his reminiscences will help younger researchers gain some perspective that will help in their studies of blazars, the most luminous long-lived objects in the universe.

The Early Years, before Jets
When the author first studied quasars in the mid-1970s, the basic theoretical model featured an expanding blob, a compact version of blast-wave scenarios borrowed from supernova remnants [3][4][5][6].Although the observed behavior agreed with the general picture of an expanding source, details of the spectral evolution differed enough from the predictions [7][8][9] to send uniform expanding spheres back to the drawing board for major modifications.One problem was that the radio spectra were too flat below the self-absorption turnover frequency, a sign of non-uniformity [10].Another discrepancy was that the peak flux of an outburst at a particular frequency did not follow the dependence on frequency of the expanding blob model [9].
The author, who was working on his Ph.D. dissertation at the time, decided to concentrate on developing and applying models involving non-uniform synchrotron sources.Unfortunately, the VLBI maps were too crude at that point to reveal the true geometry of compact extragalactic radio sources.The data were typically modeled in terms of nearly equal double components that were separating from each other at apparently superluminal speeds.The author proposed a model of a relativistic blast wave propagating through a disk-shaped wind [11,12] that was obsolete within about a year of its publication.

Relativistic Jets in Blazars
When Blandford and Rees [2] applied their radio galaxy jet model to compact radio sources in galactic nuclei, this author expressed some skepticism (found in the questions and answers after the paper in the proceedings) that a double radio structure would be as prevalent in the jet scenario as it was in the VLBI data.Within a year, this doubt proved groundless when Readhead et al. [13] used upgraded VLBI observations (including closure phases) to discover that the compact sources do, indeed, have a jet-like morphology.In the next volume of the Astrophysical Journal, Blandford and Königl [14] developed the compact jet model further, including the idea that the apparently superluminal features in the VLBI images might be shocks propagating down the jet.
These observational and theoretical breakthroughs stimulated the author to propose his own version of a relativistic jet (see Figure 1 [15]), which combined hydrodynamical acceleration and collimation, as proposed for the Galactic Center by Reynolds and McKee [16], with the downstream conical, ballistic region of the Blandford and Königl [14] model, developed further by Königl [17].The primary contribution of the author's study was to include radiative energy losses (as well as adiabatic expansion cooling, ignored in [17]).These cause gradients in the emission profile of the jet that lead to frequency-dependent time delays.The main concept is that rapid variability occurs when electrons are accelerated in small volumes or thin sections of the jet, for example, at the base of the jet or across a shock front.The highest-energy electrons, radiating at the highest frequencies, lose most of the energy within a thin layer a short distance downstream of the acceleration site.Lower-energy electrons can radiate at the lower frequencies over a larger volume.Although (ignoring the effects of absorption) the lower-frequency outburst starts at the same time as the high-frequency flare, the latter peaks and declines first, while the former builds as the larger emitting region fills with electrons.
One important prediction that the author's relativistic jet model made was that for blazars, with jets pointing almost right at us, beaming of the radiation might prevent us from detecting much emission upstream of the location where the bulk Lorentz factor of the flow reaches its asymptotic value.However, the emission near the base of the jet in less well-aligned (or slower) radio galaxy jets should be relatively bright.The less-beamed cases of M87 [18], Cyg A [19], and others, seem to follow this prediction, while the more closely aligned radio galaxy jets in 3C 111 [20] and 3C 120 [21,22] are brightest (at millimeter wavelengths) > 0.5 pc from the black hole at a feature referred to as the "mm-VLBI core."Extrapolation to blazars implies that their mm-VLBI cores are parsecs from the black hole.In M87, the feature HST-1, ~120 pc from the black hole, has properties suggesting that it might be seen as the mm-VLBI core if the jet were better aligned with our line of sight [23].Although the cm-VLBI core is probably the surface where the optical depth to synchrotron self-absorption is roughly unity [17], the linear polarization pattern of the mm-VLBI core is similar to that expected from a conical "recollimation" shock [24,25].

Combined VLBI and X-ray Observations
The landmark study of compact synchrotron sources by Jones, O'Dell, Stein, and Burbidge [26][27][28] stressed the connection between radio to optical synchrotron radiation and X-ray to γ-ray inverse Compton scattering, the latter of which they termed "synchrotron self-Compton" (SSC) emission.As both VLBI and X-ray astronomy became more developed in 1980, the author began to combine the techniques through contemporaneous observations with the VLB Network and the Einstein satellite observatory in order to use the spectral energy distributions (SEDs) to probe the physical conditions in the jets of blazars.The first focus was on the quasar NRAO140, which was one of two flaring blazars (the other was NRAO530) previously detected in X-rays by the less-sensitive, lowerresolution HEAO-1 satellite [29].The first result [30] was that, given the compactness of the radioemitting regions as measured by VLBI, relativistic beaming was needed to avoid an SSC flux that exceeded the observed value.This result was independent of the distance to the quasar, which was important because there remained a few die-hard believers in a non-cosmological interpretation of quasar redshifts being able to solve the problem of rapid variability.The observation led to the prediction of superluminal motion, which the author then measured by following changes in the structure of the jet with VLBI [31].Although this was not surprising, it was a welcome confirmation of the effects expected in relativistic jets.It also fostered the continued synergy between radio and high-energy observations.

Multi-Waveband Light Curves
Gradients in the synchrotron emission and absorption coefficients lead to frequency-dependent volumes of the emission at different frequencies, as described above.This should cause frequencydependent time delays in the peaks of flares [15].The author was heartened by reports of coincident radio and optical flares, e.g., [32,33] and, especially, Dick Miller's finding of optical flares peaking later at longer wavelengths [34,35].However, other blazar flares did not follow this pattern, which implied that blazars were unlikely to submit completely to the basic jet models.Nevertheless, these early multi-waveband observations suggested that further development of models with gradients in physical parameters would be worthwhile, as would expansion of the frequency and time coverage of future blazar monitoring programs.
The advent of sub-millimeter and far-infrared astronomy in the early 1980s unmasked this previously inaccessible region of the electromagnetic spectrum.Of particular importance was the observation of a major outburst in 3C 273 from radio to optical wavelengths by Robson et al. [36].A peculiar property of the SED of the flare was that the self-absorption turnover frequency very slightly decreased as the flux rose.This deviated from the behavior expected under any simple scenario, since whatever increases the flux should increase the opacity and, therefore, the turnover frequency as well.The author tried all of the existing variants of compact source models, finding that only bizarre combinations of parameter dependences could produce this pattern.In collaboration with Walter

Combined VLBI and X-ray Observations
The landmark study of compact synchrotron sources by Jones, O'Dell, Stein, and Burbidge [26][27][28] stressed the connection between radio to optical synchrotron radiation and X-ray to γ-ray inverse Compton scattering, the latter of which they termed "synchrotron self-Compton" (SSC) emission.As both VLBI and X-ray astronomy became more developed in 1980, the author began to combine the techniques through contemporaneous observations with the VLB Network and the Einstein satellite observatory in order to use the spectral energy distributions (SEDs) to probe the physical conditions in the jets of blazars.The first focus was on the quasar NRAO140, which was one of two flaring blazars (the other was NRAO530) previously detected in X-rays by the less-sensitive, lower-resolution HEAO-1 satellite [29].The first result [30] was that, given the compactness of the radio-emitting regions as measured by VLBI, relativistic beaming was needed to avoid an SSC flux that exceeded the observed value.This result was independent of the distance to the quasar, which was important because there remained a few die-hard believers in a non-cosmological interpretation of quasar redshifts being able to solve the problem of rapid variability.The observation led to the prediction of superluminal motion, which the author then measured by following changes in the structure of the jet with VLBI [31].Although this was not surprising, it was a welcome confirmation of the effects expected in relativistic jets.It also fostered the continued synergy between radio and high-energy observations.

Multi-Waveband Light Curves
Gradients in the synchrotron emission and absorption coefficients lead to frequency-dependent volumes of the emission at different frequencies, as described above.This should cause frequency-dependent time delays in the peaks of flares [15].The author was heartened by reports of coincident radio and optical flares, e.g., [32,33] and, especially, Dick Miller's finding of optical flares peaking later at longer wavelengths [34,35].However, other blazar flares did not follow this pattern, which implied that blazars were unlikely to submit completely to the basic jet models.Nevertheless, these early multi-waveband observations suggested that further development of models with gradients in physical parameters would be worthwhile, as would expansion of the frequency and time coverage of future blazar monitoring programs.
The advent of sub-millimeter and far-infrared astronomy in the early 1980s unmasked this previously inaccessible region of the electromagnetic spectrum.Of particular importance was the observation of a major outburst in 3C 273 from radio to optical wavelengths by Robson et al. [36].A peculiar property of the SED of the flare was that the self-absorption turnover frequency very slightly decreased as the flux rose.This deviated from the behavior expected under any simple scenario, since whatever increases the flux should increase the opacity and, therefore, the turnover frequency as well.The author tried all of the existing variants of compact source models, finding that only bizarre combinations of parameter dependences could produce this pattern.In collaboration with Walter Gear, he then turned to the shock-in-jet scenario proposed by Blandford and Königl [14], now referred to in the gamma-ray burst community as "internal shocks".

The Shock-in-Jet Model
Relativistic jet models introduced a new geometry to the interpretation of variations in blazars: Beyond the point where the bulk Lorentz factor of the flow reaches its asymptotic value, the expansion is in two dimensions rather than three.This should also be approximately the case for shocks inside the jets.In addition, electrons are expected to be accelerated mainly at the shock front, and then suffer energy losses from both radiation and adiabatic expansion.This confinement of the particle acceleration to a thin layer, along with compression of the density and component of the magnetic field that is transverse to the shock normal, promotes rapid variability.This is so because changes in the physical conditions can occur over roughly the same time-scale as the light-crossing time of the shocked region.The sudden acceleration of electrons, followed by energy losses, led to a new mathematical description of the multi-waveband behavior of blazar outbursts [37] that could be compared with the multi-waveband light curves that were beginning to be assembled.Sketches of the emission structure of the shock and of the spectral evolution are displayed in Figure 2.
Galaxies 2016, 4, 37 4 of 12 Gear, he then turned to the shock-in-jet scenario proposed by Blandford and Königl [14], now referred to in the gamma-ray burst community as "internal shocks".

The Shock-in-Jet Model
Relativistic jet models introduced a new geometry to the interpretation of variations in blazars: Beyond the point where the bulk Lorentz factor of the flow reaches its asymptotic value, the expansion is in two dimensions rather than three.This should also be approximately the case for shocks inside the jets.In addition, electrons are expected to be accelerated mainly at the shock front, and then suffer energy losses from both radiation and adiabatic expansion.This confinement of the particle acceleration to a thin layer, along with compression of the density and component of the magnetic field that is transverse to the shock normal, promotes rapid variability.This is so because changes in the physical conditions can occur over roughly the same time-scale as the light-crossing time of the shocked region.The sudden acceleration of electrons, followed by energy losses, led to a new mathematical description of the multi-waveband behavior of blazar outbursts [37] that could be compared with the multi-waveband light curves that were beginning to be assembled.Sketches of the emission structure of the shock and of the spectral evolution are displayed in Figure 2. The first phase of the outburst-which only occurs in outbursts that are initially dominated by γ-ray or X-ray emission, as sometimes occurs in quasars-is the "Compton stage".The flux rises sharply while the self-absorption turnover frequency can decline slightly, as observed in 3C 273 [36].This occurs as the inverse Compton luminosity (and therefore energy losses) declines in favor of an increase in the synchrotron flux.Although there can be complications caused by multiple scatterings if the process is mainly SSC [38], the Compton stage should occur in a manner similar to Figure 2 if the photons that are scattered originate from outside the jet or from a relatively slowly moving section of the jet.After inverse Compton losses become unimportant, the spectrum slides to lower frequencies with only a modest, if any, change in the peak flux density during the "synchrotron stage".Within this period, synchrotron losses dominate but continually decline as the shock expands while propagating down the jet, causing the magnetic field to decrease.Both the flux and turnover frequency decrease during the last, "adiabatic stage", when expansion cooling dominates.
In the same issue of the Astrophysical Journal in which the Marscher and Gear paper appeared, Hughes, Aller, and Aller [39] presented a similar model (in the adiabatic stage) that included the evolution of the linear polarization of a shock propagating down a relativistic jet.In general, shocks partially order a turbulent magnetic field by compressing the component of the field that is transverse to the shock normal.This produces polarization from a few to tens of percent, with electric vector oriented parallel to the shock normal.The ability of these two formulations to reproduce the changes in spectrum and polarization during an outburst established the shock model as the primary paradigm for interpreting blazar variability together with the appearance of superluminal knots e.g., The first phase of the outburst-which only occurs in outbursts that are initially dominated by γ-ray or X-ray emission, as sometimes occurs in quasars-is the "Compton stage".The flux rises sharply while the self-absorption turnover frequency can decline slightly, as observed in 3C 273 [36].This occurs as the inverse Compton luminosity (and therefore energy losses) declines in favor of an increase in the synchrotron flux.Although there can be complications caused by multiple scatterings if the process is mainly SSC [38], the Compton stage should occur in a manner similar to Figure 2 if the photons that are scattered originate from outside the jet or from a relatively slowly moving section of the jet.After inverse Compton losses become unimportant, the spectrum slides to lower frequencies with only a modest, if any, change in the peak flux density during the "synchrotron stage".Within this period, synchrotron losses dominate but continually decline as the shock expands while propagating down the jet, causing the magnetic field to decrease.Both the flux and turnover frequency decrease during the last, "adiabatic stage", when expansion cooling dominates.
In the same issue of the Astrophysical Journal in which the Marscher and Gear paper appeared, Hughes, Aller, and Aller [39] presented a similar model (in the adiabatic stage) that included the evolution of the linear polarization of a shock propagating down a relativistic jet.In general, shocks partially order a turbulent magnetic field by compressing the component of the field that is transverse to the shock normal.This produces polarization from a few to tens of percent, with electric vector oriented parallel to the shock normal.The ability of these two formulations to reproduce the changes in spectrum and polarization during an outburst established the shock model as the primary paradigm for interpreting blazar variability together with the appearance of superluminal knots e.g., [40,41].Association of major millimeter-wave outbursts with superluminal knots provided strong support for this picture [42].
The synchrotron stage, during which the spectrum slides toward lower frequencies while the peak flux remains roughly the same, had potential to explain the flat radio spectra of blazars.It was already known that such spectra are "conspiracies" of multiple components in the jet that emit significantly at each observed frequency [43].The other possibility, implicit in Königl's jet model [17], is that adiabatic expansion losses of the electrons as they move down the expanding jet are offset by re-conversion of flow energy to particle energy via either turbulence or shocks.

Turbulence in Blazar Jets
Besides the major outbursts and bright superluminal knots that blazars exhibit, there are many lower-amplitude variations in flux.In 1989, Miller et al. reported optical "microvariability" on a time-scale of 1.5 h in BL Lac [44].In addition, the linear polarization varies, sometimes erratically and other times more smoothly.A number of rotations of the polarization vector were reported in the early 1980s at both radio and optical wavelengths [45][46][47].It was soon recognized that the erratic variations, and even the smooth polarization rotations, could result from turbulent (or otherwise disordered) magnetic fields in the jet flow [48,49].The author, with Gear and graduate student John Travis, added this complication, as well as bending of the jet, to the shock model in an attempt to provide a more comprehensive explanation of blazar variability [50].
One question regarding turbulence is whether it merely represents the local "weather" of the plasma or plays an important role in a blazar jet.Second-order Fermi particle acceleration, as well as multiple magnetic reconnections, are a natural consequence of chaotic magnetic fields.The author maintains that these processes fill the ambient jet with relativistic particles, so that shocks only need to increase the energy per particle by about one order of magnitude in order to cause a high-energy flare.It can also cause diffusive particle acceleration when a turbulent cell with a favorable magnetic field direction (nearly parallel to the shock normal [51]) crosses a shock.This would occur erratically, causing sporadic high-energy flares such as those seen in blazars.

The "Core" of a Blazar Jet
For various reasons, there has been considerable confusion with the blazar community over the so-called "core" seen on VLBI images.Some automatically associate it with the central engine, thinking that the core is the base of the jet, in the immediate environs of the black hole.While this appears to be the case in some radio galaxies (see Section 3), there is considerable evidence that in jets that are pointing within ~15 • of the line of sight, the mm-VLBI core is >0.5 pc from the black hole.This can be inferred from the continuum spectrum, which becomes completely optically thin in the 30-300 GHz range, a fact that constrains the compactness and therefore size of the core [52].This was confirmed through combined X-ray and VLBI monitoring of the radio galaxies 3C 111 and 3C 120, with significant short-term drops in X-ray flux from the accretion disk/corona system followed about two months later by the appearance of a new superluminal knot passing through the core at 43 GHz.This travel time requires at least 0.5 pc distance between the disk and the core [20][21][22].
While Blandford and Königl [14] imagined the core to represent the point in the jet where the optical depth to synchrotron self-absorption becomes ~1, this appears to be true only at frequencies below ~30 GHz.At 43 GHz, VLBA images reveal structure in the linearly polarized intensity of cores that should not be seen if the optical thickness is the main process affecting the brightness distribution.Instead, the pattern of the polarization matches that expected from standing conical "recollimation" shocks [24,25]; see Figure 3. Since there are often one or more stationary features downstream of the core in the 22 and 43 GHz images of blazar jets [53], it is likely that the core at lower frequencies is either the location where the optical depth at that frequency ~1 or, if it is present and strong enough, the first standing shock downstream of that location [54].There are other possibilities as well: The core in a jet pointing nearly at us (i.e., in a blazar) could be the location at the end of the jet's acceleration and collimation zone [15].Alternatively, it could be where the jet bends from a larger to a smaller angle to the line of sight [55].Finally, the core could represent the most upstream location in the jet where electrons are accelerated to ultra-relativistic energies by some other means than via a standing shock.For example, the onset of turbulence owing to a current-driven instability [56] could cause acceleration by the second-order Fermi mechanism or magnetic reconnections in this region.
Galaxies 2016, 4, 37 6 of 12 upstream location in the jet where electrons are accelerated to ultra-relativistic energies by some other means than via a standing shock.For example, the onset of turbulence owing to a current-driven instability [56] could cause acceleration by the second-order Fermi mechanism or magnetic reconnections in this region.Yellow contours: Total intensity, with levels 0.5%, 1%, 2%, 4%, 8%, 16%, 32%, and 64% of the peak; color: Polarized intensity; yellow line segments denote direction of the polarization E-vector.Note the similarity between the simulated and actual images both in the pattern of the E-vectors and the upstream displacement of the peak of the polarized relative to the total intensity.(The simulated image has been rotated so that its orientation is similar to the jet of Mkn501.)

Probing Blazar Jets with Polarization and Multi-Waveband Monitoring
Magnetic fields are important to the synchrotron radiation that we see at radio to optical frequencies-and even X-ray energies in some blazars-and probably to the formation, collimation, and acceleration of relativistic jets [58,59].As Figure 3 indicates, polarization VLBI, pioneered in the mid-1980s by Bill Cotton, John Wardle, and Dave Roberts, can probe the physical conditions in compact jets [60,61].After the VLBA was commissioned in 1995, it became possible to monitor changes in the polarization structure to gain information about the dynamics of the magnetic field.Svetlana Jorstad led an international collaboration that monitored 15 blazars for three years with bimonthly VLBA images at 43 GHz, as well as measurements of the polarization integrated over each entire source at mm, submm, and optical wavelengths [62,63].There were indications that the mmsubmm-wave and optical EVPAs were often similar despite rapid variability in the optical, a result in agreement with previous observations of fewer objects by Gabuzda, Sitko, and Smith [64].Better time sampling for the quasar 0420-014 confirmed that the variable optical emission coincided with the 43 GHz core [65].
This close connection of the optical and mm-wave polarization implies that the emission regions are at least partially coincident.Most likely the optical radiation comes from a subsection of the mmwave emission region, as in the shock-in-jet model (Figure 2 [37]).The implication is broader, though: Since electrons with energies sufficient to produce optical synchrotron radiation can also scatter optical-IR photons to γ-ray energies, and since those that produce mm-wave synchrotron emission can scatter mm-IR wavelength photons to X-ray energies (see [66]), the combination of multi-epoch VLBA imaging, VLBA plus optical (and IR if attainable) polarimetric monitoring, and light curves from radio to γ-Figure 3. Left: Simulated millimeter-wave image of a conical standing "recollimation" shock, made with the author's turbulent extreme multi-zone (TEMZ) code [57], in which turbulent plasma crosses the shock, becoming compressed and energized.Black contours represent total intensity, red contours polarized intensity; black line segments are the linear polarization E-vectors.Contours are 1%, 2%, 4%, 8%, 16%, 32%, 64%, and 95% of the peak intensity.Right: VLBA image of the BL Lac object Mkn501 at 43 GHz by the author and collaborators (program VLBA-BU-BLAZARS; see web page www.bu.edu/blazars/VLBAproject.html) in July 2015.Yellow contours: Total intensity, with levels 0.5%, 1%, 2%, 4%, 8%, 16%, 32%, and 64% of the peak; color: Polarized intensity; yellow line segments denote direction of the polarization E-vector.Note the similarity between the simulated and actual images both in the pattern of the E-vectors and the upstream displacement of the peak of the polarized relative to the total intensity.(The simulated image has been rotated so that its orientation is similar to the jet of Mkn501.)

Probing Blazar Jets with Polarization and Multi-Waveband Monitoring
Magnetic fields are important to the synchrotron radiation that we see at radio to optical frequencies-and even X-ray energies in some blazars-and probably to the formation, collimation, and acceleration of relativistic jets [58,59].As Figure 3 indicates, polarization VLBI, pioneered in the mid-1980s by Bill Cotton, John Wardle, and Dave Roberts, can probe the physical conditions in compact jets [60,61].After the VLBA was commissioned in 1995, it became possible to monitor changes in the polarization structure to gain information about the dynamics of the magnetic field.Svetlana Jorstad led an international collaboration that monitored 15 blazars for three years with bimonthly VLBA images at 43 GHz, as well as measurements of the polarization integrated over each entire source at mm, submm, and optical wavelengths [62,63].There were indications that the mm-submm-wave and optical EVPAs were often similar despite rapid variability in the optical, a result in agreement with previous observations of fewer objects by Gabuzda, Sitko, and Smith [64].Better time sampling for the quasar 0420-014 confirmed that the variable optical emission coincided with the 43 GHz core [65].
This close connection of the optical and mm-wave polarization implies that the emission regions are at least partially coincident.Most likely the optical radiation comes from a subsection of the mm-wave emission region, as in the shock-in-jet model (Figure 2 [37]).The implication is broader, though: Since electrons with energies sufficient to produce optical synchrotron radiation can also scatter optical-IR photons to γ-ray energies, and since those that produce mm-wave synchrotron emission can scatter mm-IR wavelength photons to X-ray energies (see [66]), the combination of multi-epoch VLBA imaging, VLBA plus optical (and IR if attainable) polarimetric monitoring, and light curves from radio to γ-ray frequencies can allow one to map the multi-waveband emission in a blazar jet.This assumes, of course, that the blazar being monitored cooperates!However, they apparently do, at least some of the time.This is the basis of the multi-waveband monitoring program that the author and his international collaborators have pursued since shortly before the launch of Fermi and the resultant availability of continuous γ-ray light curves since August 2008.Early multi-waveband monitoring in the 1990s, when blazars could only be observed for 1-3 weeks at a time with the less sensitive EGRET detector of the Compton Gamma Ray Observatory, found that γ-ray flares tend to follow the onset of mm-wave outbursts and passage of new superluminal knots through the mm-wave core [67,68].Recent results with the Fermi data confirm this, leading to the conclusion that at least 2/3 of all γ-ray flares occur on parsec scales [69].A number of theorists have gone into denial over this result, since blazar SEDs conform with models in which jet electrons scatter optical emission-line or hot dust photons to γ-ray energies, but these seed photon fields are thought to be weak parsecs from the black hole along the jet direction.In addition, it is easier to understand how variability time-scales of flares can be as short as ~10 min if the flares occur near the black hole.There is, though, corroborating evidence: Extremely rapid very high-energy (VHE) γ-ray flares that could not escape the sub-parsec region and its high radiation field without severe attenuation from pair production e.g., [70].In addition, a mega-outburst in 3C 454.3 had a narrow simultaneous peak at mm to γ-ray frequencies, but the 1 mm emission would be optically thick if it were squeezed inside the broad emission-line region [71].Alternative proposals, such as the presence of stray emission-line [72] or dust clouds near the jet, have not been examined in detail, although emission-line variability during γ-ray flares has been reported [72,73].

Comparison of Jet Simulations with Observations
Since nature seems to have no problem making powerful relativistic, well-collimated jets with highly variable emission, we should be able to reproduce this in computer models based on fundamental physical principles.Time-variable emission from a relativistic jet was simulated by José Luis Gómez in the mid-1990s, using a relativistic hydrodynamical code developed by Ibañez and Martí [74].Moving shocks caused by velocity fluctuations at the input point created superluminal radio knots, and the changes in pressure resulted in the formation of standing conical shocks, as seen in actual jets.In addition, magneto-hydrodynamic (MHD) simulations of jets within about 10,000 gravitational radii of the black hole have provided a way of testing the magnetic launching scenario [58,59,[75][76][77][78]: A strong helical field combined with rotation of the flow might cause rotations of the polarization vector of any emission arising from this zone [79].The author's collaboration found such a rotation in the optical polarization vector of BL Lac when a new superluminal knot was detected slightly upstream of the core in 43 GHz VLBA images [80].They observed an even more impressive rotation by 720 • during a major γ-ray outburst in the quasar PKS 1510-089 [81].These observations included γ-ray flares that occurred when new superluminal knots were passing through the mm-VLBI core as well as before they reached the core.This implies that there are multiple γ-ray emission zones, some upstream and some within the core.(Other observations find γ-ray flares downstream of the core as well; e.g., [82].) The author has developed an emission code that simulates turbulent plasma crossing a standing conical shock in a blazar jet [57].The numerical model reproduces the general characteristics of blazar variability: Randomly occurring major outbursts seen across wavebands but with short-term fluctuations that often appear at optical but not γ-ray frequencies, and vice versa, plus linear polarization that fluctuates rapidly, usually about a mean value and mean position angle.The author plans to use the code to test various particle acceleration scenarios by comparing the statistical behavior of simulated versus observed variations.

Potential Solutions to Blazar Conundrums
The outstanding problems of blazars have been with us for decades: Extreme variability on time-scales shorter than we can understand, intraday flares occurring parsecs downstream of the black hole, brightness temperatures beyond the ~10 12 K inverse Compton limit [83] (measured with VLBI including the space-based RadioAstron antenna [84]), the absence of an obvious source of seed photons for producing γ-ray on parsec scales, and acceleration of particles to energies >30 GeV on time-scales of minutes, then cooling them just as rapidly.One way to solve the time-scale problem is with Doppler factors ~100, which once would have been considered beyond belief, but now seems rather conservative compared to γ-ray burst models.Nevertheless, there is a serious issue with the required high bulk Lorentz factors: The relativistic beaming is so narrow that a parent population exceeding the supply of available galaxies would be needed given the number of blazars that would require this solution, see [85].A possible way to avoid this is for a given jet to have a range of Lorentz factors over different streamlines [86].Another would be for the ultra-high Lorentz factors to occur occasionally rather than continuously.This is the solution offered by Giannios et al. [87], who have proposed that magnetic reconnections produce ultra-fast "mini-jets" of plasma with highly energetic electrons accelerated during the reconnection.This could possibly explain the strong γ-ray emission observed in the mis-aligned radio galaxy 3C 84 [88].However, the magnetic reconnection process in relativistic plasmas has yet to be studied sufficiently to determine whether mini-jets with such high Lorentz factors can be produced with sufficient coherence to explain sudden flares, or whether the travel time of magnetic realignments is too slow [89].Relativistic turbulence, with occasional ultra-fast cells might provide similar benefits as mini-jets [90], although simulations of the dynamics of such a plasma are also lacking.In support of this scenario, on the other hand, are the author's computations of variable multi-waveband flux and polarization in a many-zone model [57], which resemble observations of blazars [70,90].
In general, the parent population problem could be overcome if jets either (1) contain sporadic mini-jets in different directions that very narrowly beam their radiation over a range of directions at different times; (2) have ultra-high Lorentz factors only a small fraction of the time; or (3) always have ultra-fast spines, but these only occasionally contain enough particles (or sufficiently strong magnetic fields) to radiate significantly.More theoretical work is needed to determine which scenarios can satisfy one of these criteria as well as agree with observational data.
In the author's opinion, particle acceleration in blazar jets is a cooperative venture in which second-order Fermi acceleration and magnetic reconnections in the turbulent jet plasma provide an ambient population of relativistic particles with a power-law energy distribution.High-energy flares then occur when shocks compress and further energize these already energetic particles.Perhaps some magnetic reconnections can produce major flares as well.
The author's scientific lineage includes his advisor, Bob Brown (who co-discovered Sgr A*), his advisor Bob Gould (theorist at UCSD), well-known astrophysicist Ed Salpeter, and physicist Rudolph Peierls, whose advisors were Wolfgang Pauli and Werner Heisenberg.Of course, Heisenberg is famous for his Uncertainty Principle.Some of the intellectual aspects of this "ancestry" seems to have passed on through the generations to the author, who remains rather uncertain regarding the nature of blazar jets!

Galaxies 2016, 4 , 37 3 of 12 Figure 1 .
Figure 1.Sketch of the relativistic jet model proposed by the author in 1980 [15].Here, Γ is the bulk Lorentz factor of the flow, r is the cross-sectional radius, and R is the distance from the base of the jet.

Figure 1 .
Figure 1.Sketch of the relativistic jet model proposed by the author in 1980 [15].Here, Γ is the bulk Lorentz factor of the flow, r is the cross-sectional radius, and R is the distance from the base of the jet.

Figure 2 .
Figure 2. Left: Sketch of the Marscher and Gear shock model [37].The emitting volume increases with decreasing frequency.Right: Sketch of the evolution of the continuum spectrum according to the model.Stage 1 occurs when inverse Compton energy losses dominate, stage 2 when synchrotron losses dominate, and stage 3 when adiabatic expansion is the primary cooling mechanism.

Figure 2 .
Figure 2. Left: Sketch of the Marscher and Gear shock model [37].The emitting volume increases with decreasing frequency.Right: Sketch of the evolution of the continuum spectrum according to the model.Stage 1 occurs when inverse Compton energy losses dominate, stage 2 when synchrotron losses dominate, and stage 3 when adiabatic expansion is the primary cooling mechanism.

Figure 3 .
Figure 3. Left: Simulated millimeter-wave image of a conical standing "recollimation" shock, made with the author's turbulent extreme multi-zone (TEMZ) code [57], in which turbulent plasma crosses the shock, becoming compressed and energized.Black contours represent total intensity, red contours polarized intensity; black line segments are the linear polarization E-vectors.Contours are 1%, 2%, 4%, 8%, 16%, 32%, 64%, and 95% of the peak intensity.Right: VLBA image of the BL Lac object Mkn501 at 43 GHz by the author and collaborators (program VLBA-BU-BLAZARS; see web page www.bu.edu/blazars/VLBAproject.html) in July 2015.Yellow contours: Total intensity, with levels 0.5%, 1%, 2%, 4%, 8%, 16%, 32%, and 64% of the peak; color: Polarized intensity; yellow line segments denote direction of the polarization E-vector.Note the similarity between the simulated and actual images both in the pattern of the E-vectors and the upstream displacement of the peak of the polarized relative to the total intensity.(The simulated image has been rotated so that its orientation is similar to the jet of Mkn501.)