Why would a tidally-locked rocky planet have a first-order spherical harmonic surface temperature distribution? The new Letter to Nature Absence of a thick atmosphere on the terrestrial exoplanet LHS 3844b (also ArXiv) analyzes the thermal infrared light curve from the system (about 4.5 to 5.5 um). The planet is assumed to be tidally locked, so lack of asymmetry in the curve is cited as evidence that there is not thermal inertia due to a thick atmosphere, which is what one would expect for this planet.

In the beginning of the paper the authors say:


  We fit the extracted light curve with a simultaneous model of the astrophysical signal and the instrument behavior. The astrophysical signal consisted of a transit model and a first-degree spherical harmonics temperature map to represent the planet’s thermal phase variation.


and later:


  In addition to the spherical harmonics model, we also tested a sinusoid model, which has been commonly used to fit other phase curve data. 


I am thinking that the incident flux at a given point on the tidally-locked planet would be 

$$I \sim \max (0, \cos(\theta))$$

where $\theta$ is the static zenith angle at a given point, and so the temperature would be something like

$$T \sim I^{1/4} \sim \max (0, \cos(\theta))^{1/4}.$$

Why do they use a first-order spherical harmonic model instead? Is it related to the thermal conductivity of rock?
It probably wouldn't. But when studying these things, you don't want to go in assuming you know more than you do or you might bias the analysis.

I am thinking that the incident flux at a given point on the tidally-locked planet would be
$$I∼\max(0,\cos(\theta))$$
where θ is the static zenith angle at a given point, and so the temperature would be something like
$$T∼I^{1/4}∼\max(0,\cos(\theta))^{1/4}.$$

Unfortunately this assumes that the planet is a blackbody with zero heat redistribution over the surface, e.g. by winds or currents in a magma ocean, and no sources of heat on the nightside such as tidally-driven hyperactive volcanism which is something you certainly do not know starting out.
Spherical harmonics are a generic set of basis functions over a sphere, so it does make sense as a fit which does not assume any physical processes operating or their relative importance. In fact, this is noted in Louden &amp; Kreidberg (2018) &quot;SPIDERMAN: an open-source code to model phase curves and secondary eclipses&quot; (reference 15), which is referenced when they talk about using a spherical harmonic model. Conveniently, the paper shares the lead author with the LHS 3844 b paper, so presumably this does reflect some of the thought that went into the LHS 3844 analysis. A relevant quote from that paper:

A useful and physics-independent model is a sum of spherical harmonics. This method was used for the case of the phase curve of HD 189733b by Majeau et al. (2012). An example map generated by SPIDERMAN is displayed in Figure 4. The main observational features of a phase curve, including the offset hotspot, can typically be recovered with a only the first spherical harmonic, with the centre offset from the substellar point. (Cowan et al. 2017) explore the effects of odd harmonics in phase curve data, and find that these can correspond to weather features in the planet atmosphere.

(emphasis mine)
Kreidberg et al. (2019) do note that a simple sinusoidal fit produces unphysical negative temperatures on the nightside which can be corrected using odd harmonics, hence their choice of the first-order model.
Physics-based models probably would start to introduce a whole bunch of poorly-constrained parameters, probably overkill for the first analysis: this is mapping the terrain, figuring out how the terrain ended up like that is another stage entirely. This is a Nature paper, so brevity is par for the course.