Light curves and polarizations of gravitationally lensed kilonovae

It is anticipated that binary NS and BH-NS mergers will release neutron-rich ejecta. Here, we consider that the ejecta is composed of red and blue components with a fast velocity tail in the outermost layer of both. We assume that the half-opening angle of the lanthanide-rich red component in the equatorial region is $30^{\circ}$, and the rest of the remaining area is the lanthanide-free blue component (Bulla et al., 2019). The red or blue component, as a basic component, has a power-law density profile that can be expressed as (e.g., Hotokezaka et al., 2013; Nagakura et al., 2014; Matsumoto, 2018)

where $M_{\rm{base}}$, $\beta_{\rm{b}}=4$, $v_{\rm{b}}=0.1~{}c$, and $v_{\rm{t}}=0.3~{}c$ are the mass of red or blue component, the power-law index, the lowest velocity of the basic component, and the lowest velocity of the fast tail. We set $M_{\rm base}=0.03~{}M_{\odot}$ and $0.02~{}M_{\odot}$ for the fiducial values, respectively (for a review see Metzger, 2019).

The ejecta is close to homologous expansion after mergers, so the fast tail moving ahead in the $r$-process ejecta has the density profile (e.g., Kyutoku et al., 2014; Hotokezaka et al., 2018)

where $M_{n}=10^{-4}~{}M_{\odot}$, and $\beta_{\rm{t}}=6$ are the total tail mass (see below, for the composition of the tail), and the power-law index of the fast tail, respectively.

Since the material distribution is anisotropic, we calculate the material distribution of a unit solid angle. The mass of matter outside the $v$-shell is

Regarding the chemical composition of the different layers, we assume that the fastest layer ($v>v_{\rm n}$, $v_{\rm{n}}=0.5~{}c$ is the lowest velocity of the free neutrons) contains only free neutrons, i.e., $X_{\mathrm{r}}(v,t)=0$, and the red or blue component ($v_{\rm{b}}<v<v_{\rm{t}}$) contains only the $r$-process materials, i.e., $X_{\mathrm{r}}(v,t)=1$. For the inner part of the high-velocity tail component ($v_{\rm{t}}<v<v_{\rm{n}}$), the mass fractions of the $r$-process elements, free neutrons and protons (the products of the $\beta$-decay of neutrons) are

respectively, where the index $\alpha=10$ represents how drastically the nucleon abundance decreases in the $v<v_{\rm n}$ range, and $t_{\rm n}=900~{}\rm s$ is the $\beta$-decay timescale for neutrons.

The opacity in the ejecta arises from the electron scattering and the bound-bound transition of the $r$-process elements, so the optical depth is

where $\kappa_{\mathrm{r,red}}=10~{}\rm cm^{2}~{}g^{-1}$ and $\kappa_{\mathrm{r,blue}}=1~{}\rm cm^{2}~{}g^{-1}$ are the bound-bound opacities of the red and blue components, respectively, and $\kappa_{\mathrm{es,H}}=0.04~{}\rm cm^{2}~{}g^{-1}$ is the Thomson scattering opacity of hydrogen. In Equation (7), we neglect the contribution of electrons produced by ionization of $r$-process elements because of the low opacity. We solve for the diffusion layer velocity $v_{\rm d}$ and the photospheric velocity $v_{\rm ph}$ at each time with $\tau=c/v$ and $\tau=1$, respectively.

The bolometric luminosity is provided by the total mass of $r$-process elements and free neutrons contained outside the diffusion layer, which are given, respectively, by

and

The bolometric luminosity is approximately calculated by

where $q_{\mathrm{r}}$ and $q_{\mathrm{n}}$ are the particular heating rates of the radiative decay from $r$-process elements and free neutrons, which are given by Wanajo et al. (2014) and Kulkarni (2005), respectively, as shown below:

and

The effective temperature at the radius of the photosphere, i.e., $R_{\rm ph}=v_{\rm ph}t$, is

where $\sigma_{\mathrm{SB}}$ is the Stefan-Boltzmann constant.

The observed flux per unit solid angle is

where $h$, $k_{\rm B}$ and $D_{L}$ are the Planck constant, Boltzmann constant and luminosity distance of the source, respectively.

The observation time is

where $t$ is the emission time of the photons and $\cos\alpha=\sin\theta\cos\varphi\sin\theta_{\rm obs}\cos\varphi_{\rm obs}+\sin\theta\sin\varphi\sin\theta_{\rm obs}\sin\varphi_{\rm obs}+\cos\theta\cos\theta_{\rm obs}$ is the direction cosine between the direction of a unit solid angle at the photosphere radius and the line of sight.

If we assume that the observer is located at $\varphi_{\rm obs}=0^{\circ}$, the range of longitudes visible to an observer located in a certain $\theta$ is $[-\varphi(\theta),\varphi(\theta)]$, where $\varphi(\theta)$, as the half-day arc, is given by

If the kilonova emission is horizontal axisymmetry, the total flux of that reaching the observer at $\theta_{\rm obs}$ is written as

Finally, we convert $F_{\nu}\left(\nu,t_{\mathrm{obs}}\right)$ to a monochromatic AB magnitude by the following equation:

Extensive research has been conducted to estimate the degree of polarization in nonspherical supernovae explosions. The polarization of oblate or prolate spheroidal atmospheres was specifically computed analytically by Shapiro & Sutherland (1982) as a function of the asphericity parameter. Hoflich (1991) and Kasen et al. (2003) used a different Monte Carlo method to calculate photon propagation.

We consider that a kilonova consists of two components, red and blue, with a fast tail in the outermost layer. When the photosphere is located in the fast velocity tail, electron scattering dominates the opacity and causes a large degree of polarization. When the photosphere recedes into the $r$-process element-rich ejecta, the opacity is dominated by the bound-bound transition, which does not result in polarization, as opposed to electron scattering (e.g., Kasen et al., 2013; Tanaka & Hotokezaka, 2013). Of course, the ionization of the $r$-process elements can supply electrons, the opacity of which is much lower than that of the nucleon ejecta in our calculation (e.g., Kasen et al., 2013; Kyutoku et al., 2015). Bulla et al. (2019) found that electron scattering opacities are comparable to the bound-bound opacities in the lanthanide-free component at R band and 1.5 days after the merger in their simulations.

Brown & McLean (1977) derived an analytic expression of polarization for optically thin Thomson scattering in an axisymmetric envelope. Optically thin condition is a critical approximation for avoiding multiple scattering; and the polarization level in the optically thin case is proportional to the optical depth for electron scattering. The net polarization is as follows (e.g., Brown et al., 2000; Friend & Cassinelli, 1986):

where $\tau_{\mathrm{es}}$ is the electron scattering opacity at the photosphere, which is calculated as

including the contributions of the electrons produced by the $\beta$-decay of free neutrons and supplied by the ionization of the $r$-process elements (the first and second terms), which we ignored in Equation (7). $x$ and $A$ are the degree of ionization and the mean mass number of the $r$-process elements, respectively. Following Matsumoto (2018), we set $A=80$ and $x=1$ here.

A massive BH with $M_{l}<10^{7}~{}M_{\odot}$ can be described as a PM model (Schneider et al., 1992). The Einstein angle that corresponds to the point mass lens $M_{l}$ is defined as

where $D_{s}$, $D_{l}$ and $D_{ls}$ are the angular diameter distances to the source at redshift $z_{s}$, to the lens at redshift $z_{l}$, and between the source and the lens, respectively.

The dimensionless lens equation is

where $y=\beta/\theta_{\rm E}$, $\beta$ denotes the angular positions of the source, and $x_{\pm}=\theta_{\pm}/\theta_{\rm E}$, where “+” and “$-$” denote the parity of the image and $\theta_{\pm}$ indicates the angular positions of the image as shown in Figure 1. The corresponding two solutions of the lens equation are

The positive parity image is always the leading image, so the time delay between images $x_{-}$ and $x_{+}$ is

Their magnifications are

respectively.

The SIS model describes galaxies as gravitational lenses (Schneider et al., 1992). Its Einstein angle can be written as

where $M(\theta_{E})$ and $\sigma_{v}$ are the mass within the Einstein radius and the velocity dispersion of the lens, respectively. The corresponding lens equation is

If $y<1$, there are two solutions, i.e., $x_{\pm}=y\pm 1$. However, if $y>1$, the lens equation only has one solution: $x=y+1$. We only focus on the case of two images to highlight the role of the gravitational lensing effect. The time delay between them is

The magnifications of the images are

For a point-mass + external shear lens model, i.e., a CR model, an external component of shear from a galaxy is added to the PM lens (e.g., Chang & Refsdal, 1979, 1984; Schneider et al., 1992; Chen et al., 2021). Figure 3 in Gao et al. (2022) depicts the critical curves and transition curves in the plane of the observer at external shear strength $\gamma=0.1$. The observer can perceive four images inside the critical curves and two images beyond this range. Additionally, the later arriving image is brighter inside the range of the transition curves; otherwise, the earlier arriving image is brighter. The dimensionless lens equations are

where we set $\gamma$ equal to $0.1$. The lens equations can have up to four solutions, but they are not easy to derive. Here, we show only special solutions along the axis (i.e., $y_{1}=0$ or $y_{2}=0$).

For $y_{1}=0$, we obtain

and

For $y_{2}=0$, we obtain

and

Note that Equations (32) and (34) appear only when $\left|y_{2}\right|\leqslant 2\gamma/\sqrt{1+\gamma}$ and $\left|y_{1}\right|\leqslant 2\gamma/\sqrt{1-\gamma}$, respectively. The magnification of each image is

The time delay between the two images is

where $\overrightarrow{x_{f}}$ and $\overrightarrow{x_{s}}$ are the positions of the earlier and later images, respectively, and

The light curves of the bolometric luminosity and apparent magnitude in the R-band without the gravitational lensing effect are shown in Figure 2. The black total light curve shows that the farthest free neutron heating generates a precursor at $\sim 1~{}\rm hr$, and the peak time of ejecta heating is $\sim 60~{}\rm hrs$. As shown by the orange line, the peak apparent magnitude in the R-band of the ejecta heating is $\sim 18.6$.

In Figure 3, to highlight the bump due to the effect of gravitational lensing, we use linear coordinates to prevent confusion with the intrinsic bump of the kilonova in light curves, and the solid blue line indicates the light curve without the gravitational lensing effect. The number of images, the magnification, and the time delay between images all affect the lensed kilonova light curves.

For the PM model, three different cases of lens parameters are configured as follows: ($\beta=0.5~{}\theta_{E},M_{l}=10^{3}~{}M_{\odot}$, black dashed line), ($\beta=0.5~{}\theta_{E},M_{l}=10^{7}~{}M_{\odot}$, red solid line), and ($\beta=0.8~{}\theta_{E},M_{l}=10^{7}~{}M_{\odot}$, green solid line), from which the time delay of two images is on the order of $\sim 10^{2}~{}\rm s$ or less and the magnification of the second image is smaller than that of the first image. Because the time delay between the two images is significantly less than the ejecta heating timescale, it is impossible to distinguish the flux contribution of each image. Figure 3(a) depicts a lensed kilonova that was created by superimposing two images with successive arrivals. The magnification of each image is mainly determined by $\beta$, i.e., the source angular position relative to the lens. For the cases of $\beta=0.5~{}\theta_{\rm E}$ and $\beta=0.8~{}\theta_{\rm E}$, the peak apparent magnitudes at $\sim 60~{}\rm hrs$ are $\sim 17.8$ and $\sim 18.2$, and the increase factors are $\sim 1.04$ and $\sim 1.02$ compared to the case without the gravitational lensing effect (solid blue line), respectively. For the cases of ($\beta=0.5~{}\theta_{E},M_{l}=10^{3}~{}M_{\odot}$) and ($\beta=0.5~{}\theta_{E},M_{l}=10^{7}~{}M_{\odot}$), the light curves are exactly the same due to the same $\beta$ value and the short time delay.

For the SIS model, the lens parameters are ($\beta=0.5~{}\theta_{E},M_{l}=10^{8}~{}M_{\odot}$, black solid line), ($\beta=0.5~{}\theta_{E},M_{l}=10^{10}~{}M_{\odot}$, red solid line), and ($\beta=0.8~{}\theta_{E},M_{l}=10^{10}~{}M_{\odot}$, green solid line). The time delay would be $\sim 10^{3}-10^{5}~{}\rm s$. As seen in Figure 3(b), the ($\beta=0.5~{}\theta_{E},M_{l}=10^{8}~{}M_{\odot}$) case has a peak apparent magnitude $\sim 17.2$ at $\sim 60~{}\rm hrs$, and it is difficult to distinguish a lens signal due to the short time delay. The ($\beta=0.5~{}\theta_{E}$, $M_{l}=10^{10}~{}M_{\odot}$) case has a peak apparent magnitude $\sim 17.4$ at $\sim 54~{}\rm hrs$ and shows a lens signal, which is a plateau period with an apparent magnitude of approximately 17.5 between $\sim 30$ and $\sim 90$ hrs. The ($\beta=0.8~{}\theta_{E}$, $M_{l}=10^{10}~{}M_{\odot}$) example also exhibits a tiny lens signal, i.e., an identifiable bump with an apparent magnitude $\sim 19$ at $\sim 115~{}\rm hrs$. The reason for the appearance of the lens signal is because the lens mass is large enough to produce a time delay longer than the ejecta heating timescale.

For the CR model, we assumed that the observer’s position in the $y_{1}-y_{2}$ plane is $[0,0.15]$ and $[0,0.3]$, inside and beyond the range of the critical curve, respectively, so four or two images are available; these positions are within the transition curve range, so the later arriving lens images are brighter than the earlier arriving images. The time delay is on the order of hundreds of seconds and below, so the different images arrive almost simultaneously. In Figure 3(c), for four images (black and red lines) and two images (green and orange lines), the peak apparent magnitude at $\sim 60~{}\rm hrs$ can reach $\sim 15.2$ and $\sim 17$, and the increase factors compared to the absence of a lens effect are $\sim 1.22$ and $\sim 1.09$, respectively. The luminosity of the kilonova can be greatly magnified, but there is no bump-like lensing signal. The CR model has brighter kilonovae than the PM model due to the presence of extra shear from the galaxy.

Based on the lensed light curve results for the PM, SIS, and CR models, we discover that the time delay induced by the SIS model is the most pronounced and results in the most discernible lens signal at $\sim 100~{}\rm hrs$, and that the CR model amplifies the light curve the most. $\beta$ determines the magnification, and the lens mass mainly determines the time delay. When the lens mass is sufficiently large, the time delay is greater than the ejecta heating peak time, and we may directly witness the relatively strong gravitationally lensed kilonova signal.

The polarization is calculated as a function of electron density by Equation (19), where polarization in the polar region is negative, and polarization in the equatorial region is positive. We underscore here that the polarization at any given time is the sum of the contributions of all components.

As shown by the blue solid lines in Figure 4, if there is no gravitational lensing effect, the initial degree of polarization is almost constant with a magnitude of $\sim-0.4\%$ during the first $\sim 2~{}\rm hrs$ and then decreases smoothly. The temporal representation of the polarization can be understood as follows. The photosphere is initially in the layer of free neutrons, $\tau_{\rm es}\simeq 1$, and the total polarization is $\sim-0.4\%$. The timescale of the initial plateau can be used to effectively estimate the total mass of neutrons (e.g., Matsumoto, 2018; Li & Shen, 2019). When the photosphere recedes into the $r$-process element-rich ejecta, the main source of opacity is not electron scattering but bound-bound transition, and the recombination of electrons to the $r$-process elements causes $\tau_{\rm es}$ to drop rapidly, resulting in a rapid decrease in polarization degree. For lensed kilonovae polarization, only the number of images and time delay can affect the degrees.

For the case of the PM model in Figure 4(a), the polarization with different lens parameters varies mainly in the first 0.2 hrs because the second image is delayed by approximately 100 seconds or less from the first image. When the second image arrives, there is an abrupt rise in the polarization degree. Within the first 2 hrs, the maximum superimposed polarization is approximately $-0.75\%$. Therefore, we can judge from the first 0.2 hrs of polarization whether the kilonova is lensed by the PM lens.

For the SIS model in Figure 4(b), if $M_{l}=10^{10}~{}M_{\odot}$, a spike-like bump lasting up to 2 hrs appears at $\sim 47~{}\rm hrs$ and $\sim 69~{}\rm hrs$ with the same polarization degree $-0.45\%$ for $\beta=0.5~{}\theta_{E}$ and $\beta=0.8~{}\theta_{E}$, respectively. This is because the time delay between the two images is approximately $10^{3}-10^{5}~{}\rm s$, and the plateau phase of the second image polarization is overlaid on the falling phase of the first image polarization. For the case of ($\beta=0.5~{}\theta_{E},M_{l}=10^{8}~{}M_{\odot}$, black solid line), the degrees of polarization are $\sim-0.4\%$ between $0.1-0.6~{}\rm hrs$ and $\sim-0.75\%$ between $0.6-2~{}\rm hrs$ and then decrease smoothly. If a time delay between two images is greater than 2 hrs, a spike-like bump lasting up to 2 hrs will be produced during the late descending evolutionary stage of polarization, but the maximum polarization will be diminished. If the time delay is less than 2 hrs, a truncated plateau-like elevation will be produced in the first 2 hrs. This can be used to judge whether the kilonova is lensed by the SIS model, and the time delay between the two images can be obtained easily.

For the CR model in Figure 4(c), the polarization superimposition of four images and two images can reach $\sim-1.52\%$ and $\sim-0.75\%$ within the first 2 hrs, respectively. For the case of four images, the polarization can still be $\sim-0.26\%$ even though the smooth descending stage reaches $100~{}\rm hrs$.

By comparing the polarization under the three lens models, we find that the time delay and the number of images are the main influencing factors that cause the polarization to vary. The polarization evolution under the PM model mainly differs in the first $0.2~{}\rm hrs$. For the SIS model, a time delay greater than 2 hrs will result in a truncated bump reaching $\sim-0.75\%$ in the first 2 hrs or a spiky bump lasting up to 2 hrs in the decline phase after 2 hrs. For the CR model, the superposition of the four images makes the total polarization increase greatly. A time delay of more than 2 hrs indicates a reduced maximum polarization.

According to the polarization results of lensed kilonova, the peculiar time variability is the only chance to recognize it. For instance, the PM model shows an abrupt change in polarization at $\sim 0.2$ hrs after the merger, but obtaining polarimetry at 0.2 hrs after the merger will be extremely challenging as it requires the kilonova to be unambiguously identified right in the first few minutes after the merger and polarimetric observations to be carried out promptly on a very faint target; and the SIS model shows spike-like bumps lasting up to 2 hrs at $\sim 47$ or $\sim 69$ hrs after the merger, the duration of these features are relatively short. So the variation in polarization of lensed kilonovae is relatively difficult to be detected by the current techniques, which typically requires long exposures on faint targets as kilonovae. After accumulating enough samples, it will be possible to distinguish whether there is a lensed kilonova with polarization.