Blazars are highly variable and display complex characteristics. A key characteristic is the flux distribution or flux PDF whose shape depends upon the form of the underlying physical process driving variability. The BL Lacertae, Mrk 421 is one of the brightest and most variable blazars across the electromagnetic spectrum. It has been reported to show hints of lognormality across the spectrum from observed histograms of radio to gamma-ray fluxes. This would imply that the underlying mechanisms may not conform to the "standard" additive, multi-zone picture, but could potentially have multiplicative processes. This is investigated by testing the observed lightcurves at different wavelengths with time-series simulations. We find that the simulations reveal a more complex scenario, than a single lognormal distribution explaining the multiwavelength lightcurves of Mrk 421.

In this article, we introduce some multi-messenger studies of blazars. Main contents contain the following parts, 1. the subclasses of blazars and the observational characteristics, 2. emission variabilities from the multi-messenger observations, 3. the physical origin for the multi-messenger studies. Based on mid-IR observations collected from the $WISE$ database, we collect two sources (BZQ J1700+6830 and BZQ J1842+6809) to analyze the mid-IR variabilities. \\ (1) For BZQ J1700+6830, the period is $P\approx$ 42 hr at the four bands ($W1$, $W2$, $W3$, $W4$); for BZQ J1842+6809, the period is $P\approx$ 123 hr at the two bands ($W1$, $W2$). \\ (2) Using the periodicities, we calculate that the central black hole mass ($BH$) should be $M = (5.44\sim33.0)\times10^8M_{\odot}$ for BZQ J1700+6830 and $M = (1.13\sim6.88)\times10^9M_{\odot}$ for BZQ J1842+6809.\\ (3) For the two sources, during the flaring stage, we find anti-clockwise variation trends based on the time-based relations built by the spectral indices ($\alpha$) and flux densities ($F$).\\ (4) We calculate the Doppler factor ($\delta$) when the sources are in the dimming or brightening stages, and find that with the frequency increasing, the Doppler factor become stronger.

Relativistic shocks are one of the most plausible sites of the emission of strongly variable, polarized multi-wavelength emission from relativistic jet sources such as blazars, via diffusive shock acceleration (DSA) of relativistic particles. This paper summarizes recent results on a self-consistent coupling of diffusive shock acceleration and radiation transfer in blazar jets. We demonstrate that the observed spectral energy distributions (SEDs) of blazars strongly constrain the nature of hydromagnetic turbulence responsible for pitch-angle scattering by requiring a strongly energy-dependent pitch-angle mean free path. The prominent soft X-ray excess (``Big Blue Bump'') in the SED of the BL Lac object AO 0235+164 can be modelled as the signature of bulk Compton scattering of external radiation fields by the thermal electron population, which places additional constraints on the level of hydromagnetic turbulence. It has further been demonstrated that internal shocks propagating in a jet pervaded by a helical magnetic field naturally produce polarization-angle swings by 180$^o$, in tandem with multi-wavelength flaring activity, without requiring any helical motion paths or other asymmetric jet structures. The specific application of this model to 3C279 presents the first consistent, simultaneous modeling of snap-shot SEDs, multi-wavelength light curves and time-dependent polarization signatures of a blazar during a polarization-angle (PA) rotation. This model has recently been generalized to a lepto-hadronic model, in which the high-energy emission is dominated by proton synchrotron radiation. It is shown that in this case, the high-energy (X-ray and $\gamma$-ray) polarization signatures are expected to be significantly more stable (not showing PA rotations) than the low-energy (electron-synchrotron) signatures.

This article shortly introduces Gaussian processes as a new approach for modelling time series in the field of blazar physics. In the second part of the paper, recent results from an application of GP modelling to the multi-wavelength light curves of the blazar PKS 1502+106 are discussed.